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ELECTRONIC  MULTIPOLE  MOMENTS  AND  POLARIZABILITIES  AND 
THEIR  APPLICATION  TO  WEAK  INTERMOLECULAR  INTERACTIONS 


By 

Tadeusz  Stanislaw  Pluta 
December  1990 


Chairman:  Rodney  J.  Bartlett 
Major  Department:  Chemistry 

Electric  field  properties  of  several  small  molecules  and  their  weak  interactions  are 
studied  using  Coupled-Cluster  and  Many-Body  Perturbation  Theory  (CC/MBPT)  method. 
After  investigating  the  behavior  of  the  Ne2  pair  polarizability  as  a function  of  R,  using 
finite-field  method  CC/MBPT,  an  improved  Coupled-Cluster  density  matrix  method  is 
used  to  calculate  electric  dipole,  quadrupole  and  octopole  moments  for  the  HF  and  H2O 
molecules.  This  permits  all  multipole  moments  to  be  computed  as  expectation  values 
rather  than  using  standard  finite-field  methods.  Moments  perturbed  by  the  presence  of  an 
external  charge  are  used  to  determine  the  components  of  the  electric  tensors  of  second- 
and  higher-rank  by  differentiating  of  the  induced  moments.  Similarly,  an  external  point 
charge  finite-field  procedure  is  employed  to  determine  the  electric  tensors  of  spherical 
systems  like  Ne  and  the  fluoride  anion  F"~  which  do  not  posses  any  permanent  electric 
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moments.  The  method  is  used  to  calculate  the  higher-order  electric  properties  at  a highly 
accurate  correlated  level  of  theory  for  atoms  and  small  molecules.  Many  of  the  higher 
order  tensors  are  determined  for  the  first  time.  In  particular,  the  F-  anion  is  shown  to  have 
an  unusually  large  second  hyperpolarizability.  When  experimental  data  are  available  the 
overall  agreement  is  favorable. 

The  higher-rank  and/or  higher-order  electric  multipoles  and  polarizabilities  deter- 
mined are  successfully  used  to  study  the  long-range  intermolecular  interactions  for 
hydrogen-bonded  systems  including  both  electrostatic  and  induction  contributions.  Ap- 
plications to  the  H2O  — F~,  H2O  — HF  , and  FHF““  systems  are  presented.  In  order 
to  extend  the  applicability  of  the  multipole  analysis  to  shorter  distances  the  Cumulative 
Atomic  Multipole  Moments  (CAMM)  method  is  applied. 


CHAPTER  I 


OVERVIEW 

Theoretical  investigations  of  atomic  and  molecular  systems  mainly  involve  deter- 
mining an  approximate  solution  of  the  Schrpdinger  equation.  Modem  quantum  chemical 
methods  for  including  a large  portion  of  electron  correlation  effects  are  capable  of  yielding 
very  accurate  answers.  However,  the  total  energy  is  not  a particulary  useful  characteristic 
of  a given  system.  Other  than  the  energy,  properties  like  the  dipole  moment  and  dipole 
polarizability  are  often  needed.  Molecular  electric  moments  and  polarizabilities,  besides 
being  of  pure  academic  interest,  are  sensitive  discriminants  of  the  quality  of  molecular 
wavefunction.  They  are  also  important  in  connection  with  intermolecular  interactions. 

It  is  well  known  that  in  the  long-range  limit  the  classical  multipole  expansion  offers 
both  a valuable  physical  understanding  of  the  interaction  process  and  a very  attractive 
computational  method.  Recently  it  has  been  shown  that  simple  electrostatic  models 
can  accurately  describe  the  geometry  and  energetics  of  hydrogen-bonded  complexes, 
yet  little  attention  has  been  paid  to  the  inclusion  of  the  induction  component.  The 
objective  of  this  work  is  to  (1)  present  a method  for  the  prediction  of  higher  moments 
and  polarizabilities  of  molecules,  (2)  apply  these  methods  for  a number  of  small  molecules 
including  the  anomalous  Ne2  polarizability  curve,  and  (3)  use  the  predicted  values 
to  construct  electrostatic  models  for  weak  intermolecular  interactions  augmented  by 
induction  effects.  From  the  point  of  view  of  quantum  chemistry  electric  tensors  are 
either  first-  or  higher-order  molecular  properties.  First-order  properties  can  usually  be 
obtained  from  expectation  value  expressions.  Recently  a density  matrix  approach  has  been 
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introduced  making  it  possible  to  calculate  the  first-order  properties  from  highly  correlated 
wavefunctions.  Since  the  majority  of  current  correlated  methods  lead  to  wavefunctions 
which  do  not  satisfy  the  Hellmann-Feynman  theorem  a straightforward  application  of  an 
expectation-type  formulae  introduces  an  error  which  needs  to  be  assessed.  The  problems 
associated  with  the  use  of  the  Coupled  Cluster  wavefunction  to  calculate  the  first-order 
properties  will  be  analyzed. 

Second-  and  higher-order  electric  properties  are  obtained  from  perturbation  theory. 
One  can  employ  either  analytical  methods  or  use  simpler  numerical  differentiation 
schemes.  For  properties  like  the  second  dipole  hyperpolarizability  7 the  numerical 
approach  becomes  difficult  to  implement  and  unreliable  due  to  loss  of  precision.  The 
analytic  approach,  besides  being  more  elegant,  is  not  plagued  by  numerical  precision 
problems.  However,  both  the  theory  and  computer  implementation  is  difficult  From  a 
purely  practical  point  of  view  it  would  be  of  great  value  if  the  order  of  differentation  in  the 
numerical  scheme  could  be  reduced.  Such  a possibility  is  considered  where  the  charge 
perturbation  method  is  used  together  with  the  Coupled  Cluster  density  matrix  method 
for  the  electric  multipole  moments.  Instead  of  using  numerically  unstable  differentiation 
with  respect  to  the  energy  one  can  use  numerical  differentiation  with  respect  to  first-order 
electric  properties.  Another  advantage  of  the  method  is  that  now  many  interesting  higher 
rank  electric  tensors,  e.g.,  octopole-dipole  polarizability,  can  be  determined.  In  order 
to  assess  the  accuracy  of  the  method  it  was  applied  to  several  simple  molecules:  HF, 
H2O,  F-  and  Ne.  Many  results  presented  have  been  calculated  at  the  correlated  level 
for  the  first  time. 

Electric  multipoles  and  polarizabilities  are  then  used  to  construct  intermolecular  po- 
tentials. The  availabilty  of  molecular  polarizabilities  allows  the  introduction  of  induction 
forces  to  purely  electrostatic  models.  The  induction  interactions  prove  to  be  much  weaker 
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than  the  electrostatic  interactions.  Although  less  important  at  equilibrium  distances,  it  is 
shown  that  induction  is  important  at  geometries  different  from  equilibrium. 

Electric  properties  of  an  individual  molecule  can  be  affected  by  intermolecular 
interactions  in  the  same  way  as  the  energy.  Chapter  IV  contains  results  of  accurate 
MBPT/CCSD  + T(CCSD)  calculations  on  the  dipole  pair  polarizability  as  a function  of 
the  interatomic  distance  for  the  interesting  Ne2  system.  Since  the  calculations  are  of 
the  supermolecule  (SM)  type  the  basis  set  superposition  error  contaminates  the  results. 
There  is  a controversy  on  how  to  eliminate  or  reduce  this  error  for  interaction  energy 
calculations.  Still  less  is  known  about  the  possible  effect  of  the  basis  set  superposition 
error  on  electric  properties  of  interacting  systems.  The  present  results  indicate  that  the  Ne2 
system  behaves  differently  from  other  rare  gas  pairs.  The  detailed  analysis  of  the  electron 
correlation  contribution  and  the  basis  set  superposition  may  be  valuable  in  constructing 
the  dipole  polarizability  functions  for  other  rare  gas  systems. 


CHAPTER  II 


REVIEW  OF  THE  THEORY  OF  INTERMOLECULAR  INTERACTIONS 

Supermolecule  Approach 

The  treatment  of  interacting  molecules  A and  B within  the  framework  of  the  su- 
permolecule (SM)  method  is  in  principle  no  different  from  the  treatment  of  the  isolated 
molecules.  In  order  to  determine  the  intermolecular  interaction  energy  AE  one  has  to 
calculate  the  difference 

AE ab  = EAb  ~ (Ea  + Efl)  (1) 

where  Eab  is  the  total  energy  of  the  system  (supermolecule)  AB  and  Ea  and  Eb  denote 
the  energies  of  the  isolated  molecules.  The  supermolecule  approach  has  gained  much 
popularity  recendy  due  primarily  to  the  availability  of  general  purpose  ab  initio  computer 
codes  which  are  able  to  provide  accurate  values  for  all  the  energies  occuring  in  (1)  for 
a variety  of  systems.  The  second  major  advantage  of  the  SM  approach  lies  in  the  fact 
that  all  components  to  the  interaction  energy:  electrostatic,  induction  and  exchange- 
repulsion  are  treated  on  equal  footing  for  all  intermolecular  distances  R.  This  differs 
sharply  from  the  perturbation  theory  treatment  (to  be  discussed  in  the  next  section) 
where  different  components  demand  different  theoretical  descriptions.  If  the  energies 
in  (1)  are  determined  by  methods  which  include  electron  correlation  the  dispersion  part 
of  the  interaction  is  also  taken  into  account. 

Sometimes  the  term  “ variation”  or  variational  method  is  used  in  the  literature  in 
reference  to  the  SM  approach  (e.g.  Hobza  and  Zahradnik,  1980).  However,  even  when 
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variational  methods  are  employed  to  calculate  all  energies  in  Eq.  (1)  (e.g.,  SCF,  Cl, 
MCSCF  etc.)  the  difference  AE  is  not  variationally  bound.  In  this  work  only  the  term 
supermolecule  will  be  used  to  charcterize  methods  attempting  to  determine  AE  directly 
from  the  definition  (1). 

For  treatment  of  weak  interactions,  the  SM  method  has  some  severe  drawbacks.  The 
first  is  due  to  the  fact  that  the  SM  calculations  yield  a net  value  of  the  interaction 
energy  so  a small  energy  is  obtained  from  the  difference  of  two  large  numbers.  Also, 
the  physical  insight  into  the  nature  of  the  interaction  process  is  lost.  Other  shortcomings 
are  related  to  finite  basis  sets  used  in  calculations.  The  best  known  effect  is  the  so-called 
basis  set  superposition  error  (BSSE).  The  electronic  structure  of  the  molecule  A (B)  is 
better  described  in  the  supermolecule  AB  basis  set,  A U B , than  it  is  when  the  isolated 
molecule  A (B)  is  calculated  with  its  own  basis  set  A (B).  This  error  is  very  hard  to 
correct  in  current  basis  set  method  when  high  accuracy  is  required. 

As  a remedy  for  the  lack  of  physical  interpretation  of  the  SM  results  several  interac- 
tion energy  partitioning  schemes  have  been  introduced.  It  is  sufficient  here  to  consider 
the  basic  ideas  of  the  most  succesful  of  these  schemes,  namely  the  Morokuma  decom- 
position (Kitaura  and  Morokuma,  1976,  Morokuma  and  Kitaura,  1980).  In  this  method 
the  supermolecule  interaction  energy  AE  , obtained  at  the  SCF  level,  is  decomposed 
arbitrarily  into  a series  of  terms 

AE  = AEetst  + AE  poi  + AEexc/l  + AEcr  + AEmjj;  (2) 

Here  the  electrostatic  interaction  energy  AEeist  is  the  energy  of  interaction  between  the 
undistorted  charge  distribution  of  molecule  A and  that  of  molecule  B.  If  by  Eo  one 
denotes  the  sum  of  the  SCF  energies  for  A and  B one  can  calculate  AEelst  as  the 
difference  AEe/s<  = Ei  - Eq  where  Ei  is  the  energy  corresponding  to  the  wavefunction 
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= <2>^  • $5  • The  symbol  'I'i  is  the  Hartree  product  of  two  undistorted  wavefunctions 
for  A and  B.  There  is  no  exchange  of  electrons  nor  polarization  of  one  molecule  due 
to  the  presence  of  the  other.  The  next  step  is  to  include  the  mutual  polarization  of  the 
electron  distribution  caused  by  the  presence  of  the  partner  molecule.  Again,  no  electron 
exchange  is  allowed  and  the  corresponding  wavefunction  #2  takes  the  following  form 
\&2  = $a  • where  $a  is  the  perturbed  SCF  wavefunction  for  A in  the  presence 
of  the  charge  distribution  B and  <J>b  corresponds  to  the  perturbed  molecule  B.  The  total 
polarization  energy  can  be  expressed  as  the  difference  AE poi  = E2  — Ei  . This  difference 
can  be  further  decomposed  into  the  sum  of  two  polarization  energies  for  molecule  A and 
B.  The  exchange-repulsion  term  is  the  energy  associated  to  the  wavefunction  'I' 3 satisfying 
the  Pauli  principle  \&3  = A($()a  ■ $°D)  . The  operator  A exchanges  electrons  between 
molecules  A and  B.  The  exchange-repulsion  part  of  the  interaction  energy  is  given  as 
the  difference  between  the  energy  corresponding  to  the  ansatz  ^3  and  Ei.  In  a similar 
spirit  one  can  construct  the  Hartree  product  of  the  wavefunction  for  A (or  B)  modified 
to  account  for  the  electron  transfer  from  A(B)  to  B(A)  'I'4  /4-*b  = &A-+B  • $5 

The  difference  between  the  SCF  interaction  energy  and  the  sum  of  all  terms  described 
above  is  the  coupling  or  mixing  term  AEmix.  This  term  is  a measure  of  the  separability 
of  the  energy  components  and  becomes  large  as  the  interactions  become  stronger. 

Using  SCF  molecular  orbitals  for  monomers  A and  B one  can  construct  the  super- 
molecule AB  Fock  matrix.  Now,  having  defined  the  blocks  in  the  supermolecule  AB 
matrix  one  solves  a set  of  modified  Hartee-  Fock  equations  by  diagonalizing  only  the 
selected  parts  of  the  Fock  matrix  with  the  rest  of  it  set  to  zero.  In  general,  the  pos- 
sibility of  decomposing  the  interaction  energy  into  physically  well-defined  terms  is  of 
great  value  both  for  the  purpose  of  better  understanding  the  nature  of  the  interaction  and 
for  more  practical  reasons.  The  flexibility  of  the  method  allows  one  to  modify  it  easily 
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e.g.,  the  change  of  monomer  geometry  due  to  interaction  can  be  readily  incorporated. 
The  Morokuma  method  and  its  various  expansion  and  modifications  have  been  employed 
frequently  to  analyze  the  energetics  of  hydrogen  bonding,  protonated  complexes  (e.g., 
Umeyama  and  Morokuma,  1976)  , bonding  in  transition  metal  complexes  (Morokuma 
and  Kitaura,  1981)  etc.  Like  all  SM  method  the  Morokuma  decomposition  is  affected 
by  the  basis  set  superposition  error. 

The  limitations  of  the  energy  decomposition  scheme  are  numerous.  First  of  all,  the 
decomposition  is  not  unique.  There  is  no  criterion  except  for  rather  vague  “ physical 
interpretation”  for  adding  new  terms.  It  was  demonstrated  that  the  energy  components  are 
more  sensitive  to  the  basis  set  than  the  total  interaction  energy.  The  energy  decomposition 
often  fails  when  applied  to  strong  interactions  so  new  terms  are  introduced  to  analyze  the 
AEmix  term.  In  conclusion,  the  decomposition  methods  can  offer  physical  interpretation 
of  considerable  help.  One  should  be  aware,  however,  of  the  arbitrariness  of  the  method 
and  its  inherent  semi  quantitative  character. 

In  the  standard  SM  calculations  the  energy  of  the  supermolecule  AB  is  evaluated 
with  the  union  of  the  two  relevant  basis  sets 

Xab  = XaUxb  (3) 

The  interaction  energy  calculated  as  the  difference  of  the  supermolecule  energy  and  the 
monomer  energies  of  A and  B can  be  written  explicitly  as 

AEab(.R)  = EAb(xa  U xb)  - E^(x>4)  - E b(xb)  (4) 

This  energy  is  contaminated  by  a systematic  error  known  as  the  basis  set  superposition 
error  (BSSE).  The  error  results  from  using  the  truncated  basis  set  expansion  for  the  system. 
In  the  limit  of  a complete  basis  set  the  BSSE  should  vanish.  In  the  supermolecule  AB 
orbitals  on  B (A)  become  available  for  A (B)  thus  leading  to  the  improvement  of  energy 
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of  A (B)  in  the  supermolecule  with  respect  to  the  energy  of  the  isolated  molecules  A 
and  B.  This  is  a non-physical  effect  and  its  occurence  and  magnitude  can  change  even 
a qualitative  description  of  an  interaction,  e.g.,  artificial  minima  can  appear  on  the  SCF 
potential  energy  curve  for  the  noble  gas  dimers  as  will  be  shown  in  Chapter  IV.  It  was 
argued  (Boys  and  Bemardi,  1970)  that  in  order  to  eliminate  the  BSSE  one  should  use 
the  so  called  counterpoise  (CP)  correction  technique.  In  this  method  the  energies  of  the 
molecules  A and  B are  calculated  with  the  full  supermolecule  basis  set  The  energy  of 
the  molecule  A (B)  calculated  with  the  basis  xa^XB  is  the  proper  reference  energy  to  be 
used  in  the  SM  computations.  The  corrected  form  of  the  basic  SM  equation  reads  now 

AE ab(R)  = EAb(xa  U xb)  ~ Ea(xa  U xb ) - E b(xa  U xb)  (5) 

The  calculation  of  the  energy  of  A (B)  is  now  identical  to  the  SM  calculation  except  for 
the  different  number  of  electrons  (only  electrons  of  A present)  and  the  nuclear  charges 
of  B (or  A)  (“ghost  orbitals")  which  are  set  to  zero.  The  energy  of  A (B)  depends  on 
the  geometry  R and  it  leads  to  longer  calculations.  One  may  notice  that  calculations  can 
be  organized  in  such  a way  that  all  two-electron  integrals  in  the  AO  basis  are  calculated 
only  once  for  a given  geometry  and  then  stored.  The  CP  correction  can  be  interpreted 
as  the  sum  of  corrections  to  monomer  energies  which  change  the  reference  energy  of  the 
molecule  A(B)  under  the  interaction 

AE  2£(JQ  = AEab(R)  + *>7  + 8CbP  = 

(6) 

AEab(R)  + {E^X/O  - EjA(xa  U Xb)}  + (Ea(xb)  - E b(xa  U xb)} 

The  CP  method  utilizing  the  full  set  of  ghost  orbitals  has  met  with  much  criticism.  This 
criticism  is  based  on  the  intuitively  obvious  violation  of  the  Pauli  principle:  the  occupied 
ghost  orbitals  on  B should  not  be  used  again  as  potential  basis  functions  for  A.  The  second 
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argument  against  using  the  full  CP  scheme  is  that  this  scheme  overestimates  the  BSSE 
leading  to  an  approximate  upper  bound  for  the  interaction  energy  (Johansson  et  al.,  1973). 

Daudey  et  al.,  (Daudey  et  al.,  1974)  were  the  first  to  introduce  in  a rigorous  manner 
the  virtual-only  CP  technique.  The  monomer  modified  CP  corrections  can  be  written 
as  follows 


6C/~V  = ea(xa)  - ea(Xa  u xt') 


(7) 


6%p  v = Eb(xb)  ~ E b(xb  U xVBrt) 

where  the  virtual  orbitals  on  B (A)  must  be  obtained  by  some  projections  from  the  full 
Xa  (B)  set.  The  projection  are  not  uniquely  defined  in  the  virtual-CP  method  thus  leading 
to  many  alternative  methods  (e.g.,  Hayes  et  al.,  1984).  The  supermolecule  calculations 
can  be  either  performed  with  the  union  basis  like  in  the  full  CP  method  or,  the  modified 
supermolecule  basis  set  can  be  employed 

XCAB~V  =XAU  Xb"  U xt‘  u XB  (8) 

The  latter  option  requires  orthogonalization  of  monomer  orbitals.  It  was  shown  that  in 
practical  calculations  the  basis  sets  (3)  and  (8)  yield  similar  results  (Van  Lenthe  et  al., 
1984).  Important  criterion  of  the  correctness  of  any  procedure  attempting  to  eliminate 
the  BSSE  is  that  any  corrections  should  vanish  in  the  limit  of  a complete  basis  sets.  The 
full  CP  method  satisfies  this  criterion  but  some  of  the  proposed  schemes  like  the  common 
virtual  space  for  A,B  and  the  supermolecule  AB  (Spiegelmann  and  Malrieu,  1980)  do  not. 
Since  rigorous  virtual-only  CP  corrections  lead  into  more  complicated  calculations  and 
they  are  not  free  from  arbitrariness  in  chosing  orthogonalization,  localization  etc.  there 
have  been  proposals  to  use  approximate  virtual-only  CP  schemes  (Fowler  and  Madden, 
1983,  Olivares  del  Valle  et  al.,  1986,  Schwenke  and  Truhlar,  1985).  The  approximate 
polarization-CP  method  consists  of  using  only  polarization  atomic  orbitals  (sometimes 
with  diffuse  functions)  as  the  ghost  orbitals  on  the  partner  molecule.  This  approximate 
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scheme  was  used,  e.g.,  in  the  calculations  of  transition  metal  clusters  (Miyoshi  et  al., 
1983)  where  the  size  of  the  basis  set  is  of  paramount  importance.  Extensive  SCF  study 
of  the  (HF)2  dimer  (Schwenke  and  Truhlar,  1985)  shows  that  the  BSSE  as  evaluated  by 
the  CP  method  does  not  decrease  systematically  with  the  size  and  quality  of  the  monomer 
basis  set  so  the  resulting  energy  does  not  “converge"  to  the  correct  (but  still  unknown  !) 
value.  Application  of  the  polarization-CP  technique  saves  computational  time  but  does 
not  remove  irregularities  from  the  basis  set-BSSE  relation.  The  results  of  this  study  were 
used  as  an  argument  against  using  any  form  of  BSSE  correction.  Additional  justification 
for  this  radical  treatment  may  come  from  the  approximate  cancellation  of  the  BSSE  and 
the  part  of  the  dispersion  energy  when  polarization  functions  are  missing  from  the  basis 
sets  (Bphm  and  Ahlrichs,  1985).  Such  a cancellation  of  errors  may  occur  if  the  BSSE 
depends  on  R as  the  long-  range  (missing)  part  of  the  dispersion  energy,  i.e.  R- 6,  R- 
etc.  It  is  of  course  not  to  be  expected  for  many  systems. 

The  BSSE  cannot  be  expected  to  be  small  when  correlation  methods  are  applied 
either.  To  the  contrary,  evidence  shows  that  for  most  cases  the  correlated  part  of  the 
BSSE  may  be  greater  and  at  the  same  time  more  difficult  to  handle.  Unlike  for  isolated 
molecules  (monomers)  the  full  Cl  calculations  offer  no  definitive  answers  for  the  full  Cl 
calculations  also  contain  BSSE. 

Despite  many  papers  published  on  the  controversy,  the  question  of  how  to  eliminate 
the  BSSE  from  the  SM  results  seems  far  from  being  resolved.  The  full  CP  method  is 
most  staightforward  and  nonambiguous,  but  the  possible  violation  of  the  Pauli  principle 
cannot  be  ignored.  On  the  other  hand  the  virtual-only  CP  methods  are  arbitrary  and 
lead  to  more  complicated,  longer  calculations.  Polarization-CP  methods  cannot  be  used 
but  as  only  practical  solutions,  since  there  is  no  justification  for  them.  Calculations 
employing  no  BSSE  correction  at  all  often  yield  good  agreement  with  the  experimental 
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data.  However,  the  agreement  is  fortuitous  as  it  depends  heavily  on  the  basis  set  functions 
used  in  calculations  (Diercksen  et  al.,  1975). 

It  should  also  be  noted  that  Eq.  (1)  holds  only  for  size-extensive  methods  like 
Coupled-Cluster  theory  and  Many-Body  Perturbation  Theory.  Any  method  aiming  at 
the  elimination  of  the  basis  set  superposition  effects  must  not  destroy  this  property.  The 
BSSE  corrections  for  the  CISD  method  taking  into  account  its  non  size-extensivness  have 
been  derived  (Price  and  Stone,  1979)  but  never  have  gained  any  popularity. 

Finally,  all  the  above  described  methods  of  dealing  with  the  BSSE  can  be  classified 
as  a posteriori  adjustments,  i.e.  first  we  introduce  the  error  then  we  try  to  eliminate  or 
reduce  it.  Mayer  (Mayer,  1983,  Mayer  and  Suijan,  1989a,  1989b)  advocated  a different 
approach  which  can  be  called  a priori  elimination  of  the  BSSE.  In  order  to  do  that 
he  introduced  a non-  hermitian  Fock  operator  for  the  supermolecule  AB  with  explicit 
partitioning  into  A and  B intramolecular  parts  and  an  intermolecular  interaction  term. 
The  monomer  orbitals  were  chosen  to  be  SCF  orbitals  independently  generated  for  A 
and  B,  and  thus  nonorthogonal.  The  SCF  results  are  close  to  the  full  CP  scheme,  but  as 
one  expects,  the  correction  is  now  lower.  The  second  major  group  of  methods  directly 
avoiding  the  BSSE  is  Valence  Bond  theory  (Wormer  et  al.,  1975,  Raimondi,  1984).  Basis 
set  superposition  error  does  not  only  change  the  interaction  energy,  it  also  affects  other 
molecular  (monomer)  properties.  The  basis  superposition  effect  on  electric  molecular 
properties  will  be  discussed  in  Chapter  IV. 

The  error  introduced  by  the  basis  superposition  is  definitely  the  most  serious  one  and 
the  most  difficult  to  remove.  However,  other  errors  in  the  SM  interaction  energy  can  occur 
due  to  the  approximate  nature  of  the  basis  sets  used  in  the  calculations.  The  problem  of 
choosing  or  generating  the  appropriate  basis  functions  for  interaction  energy  calculations 
have  been  adddressed  in  the  literature.  Without  going  into  details,  one  can  state  that 
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the  requirements  that  the  basis  sets  must  meet  depend  on  the  nature  of  the  interaction 
energy.  Long-range  parts  of  the  interaction  energy  demand  an  accurate  representation 
of  electric  multipole  moments  and  multipole  polarizabilities,  while  short-range  forces 
require  an  adequate  description  of  the  valence  electron  density  region.  In  very  accurate 
studies  on  small  Van  der  Waals  systems  additional  basis  functions  have  been  located 
in  the  center  of  bonds.  The  benchmark  calculations  on  He2  (Van  Lenthe  et  al.,  1988) 
employed  up  to  130  basis  functions  with  optimized  exponents  located  on  and  off-nuclei. 
Success  of  the  model  potential  SM  calculations  (Andzelm  et  al.,  1985)  suggests  that  the 
core  electrons  are  not  of  utmost  importance.  Results  of  Schwenke  and  Truhlar  and  many 
other  studies  indicate  strongly  that  the  total  monomer  energy  should  not  be  used  as  a 
criterion  in  selecting  basis  functions. 


Perturbation  Theory  Of  Intermolecular  Interactions 

Since  the  intermolecular  interaction  energy  is  typically  of  the  order  10“ 3 to  10-8 
times  the  total  energy  of  the  system  it  seems  quite  natural  to  choose  perturbation  theory  as 
a method  for  its  theoretical  description.  This  choice  however,  leads  to  serious  problems 
at  short  distances. 


Let  us  consider  two  interacting  systems  (molecules)  A and  B with  nA  and  ns  electrons, 
respectively.  We  assume  that  the  solutions  of  the  Schrpdinger  equation  for  individual 
molecules  are  known 


*2*1  = E M 
H°b*b  = E 


(9) 


The  interaction  energy  operator  (perturbation)  V is  the  sum  of  Coulombic  interactions 
between  pairs  of  charged  particles  within  molecules  A and  B.  Explicitly,  the  form  of  the 
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operator  V is  given  as  follows 


Na  nB 


Nb  nA 


A "■B  ry  11 1!  '"A  ry  nA  -%  *•  A -■ a 

1'  = -EE--EE-+EE1+EE 

Z— / Z j r . Z — j Z — > r rj  Z— / Z — / r . • Z — y Z — ✓ 

^ 1 •_  i D _ 1 •_  1 iB  •_  i *i  /i 1 D 1 


7Vb 


>1=1  J=1  r->A 


5=1  i=l 


.=1  j=l  ' ,J 


A=1  5=1 


ZaZb 

rAB 


(10) 


The  first  two  sums  represent  intermolecular  electron-nuclear  attraction,  the  third  term  is 
a genuine  intermolecular  two-electron  repulsion  and  the  last  term  is  a constant  (for  a 
given  geometrical  arrangement  of  an  interacting  system)  nuclear-nuclear  repulsion  term. 
It  should  be  emphasized  that  all  terms  in  the  above  expression  are  of  inter  molecular 
character.  Our  goal  is  to  find  the  change  in  the  energy  of  the  system  when  molecules 
interact,  or  in  other  words  when  the  operator  V is  switched  on.  Without  V the  Hamiltonian 
for  the  system  is  just  the  sum  of  the  individual  Hamiltonians  for  A and  B 

H0  = H°A  + Hng  (11) 


The  Schrpdinger  equation  for  the  non-interacting  system  becomes  very  simple 


H°\ J>°  = E°\l>0 


(12) 


where  is  the  product  of  two  unperturbed  wave  functions  for  A and  B,  and  E°  is  the 


sum  of  unperturbed  energies  of  A and  B 


S'0  = 


E°  = + E°b 


(13) 


The  Schrpdinger  equation  for  the  interacting  systems  can  be  written  as 


(ifu  + F)^°  = E°^° 


(14) 


Having  defined  Ho,  and  V one  can  follow  the  standard  technique  of  Rayleigh- 
Schrpdinger  perturbation  theory. 

The  first-order  correction  to  the  energy  is  expressed  as  the  expectation  value  of  the 
operator  V 


Eli  = 


(15) 
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In  this  formula  an  obvious  generalization  to  include  the  interaction  between  the  states 
other  than  the  ground  state  was  introduced. 

In  second-order  perturbation  theory  one  obtains  the  following,  well-known  formula 


e'2>  = - E 


- „ E"  - E°,  + ES  - E° 

n,m  j^O  A A B B 


(16) 


(17) 


The  situation  when  both  n and  m indices  simultaneously  take  the  value  0 is  excluded  from 
this  summation.  But  nothing  prevents  us  from  restricting  one  of  the  indices  to  be  zero 
and  summing  over  the  second  one.  The  relevant  formula  consists  of  two  parts  and  reads 

bS— E 

m 

This  corresponds  to  the  induction  part  of  the  second-order  interaction  energy.  The  second 
part  of  the  general  second-order  energy  is  obtained  when  neither  the  n nor  m index  is 
allowed  to  be  zero.  This  situation  corresponds  to  the  dispersion  energy 


gm  gO  / v jTjra  g0 

B B n A A 


e(2)  = _ V*  v A 

disp  / j tpn 


^;E5-E» +es-e» 


(18) 


The  dispersion  energy  corresponds  to  the  interaction  between  two  mutually  induced 
electron  distributions.  It  has  no  classical  analogue,  and  for  the  ground  state  it  is  always 
negative. 

The  crucial  assumption  in  the  above  derivation  is  that  of  the  product  form  of  the 
wavefunction,  i.e.  'D°  = . This  particular  form  of  the  wavefunction  forces 

some  electrons  to  “belong”  to  the  molecule  A,  while  the  rest  of  them  “belong”  to  B.  The 
product  form  of  the  wavefunction  ignores  the  electron  exchange  between  A and  B.  It  is 
known  that  the  Coulombic  forces  are  long-range  as  opposite  to  the  short-range  exchange 
interaction.  This  difference  in  range  makes  this  version  of  perturbation  theory  still 
applicable  at  the  large  intermolecular  distances.  The  name  Polarization  Approximation 
(PA  , Hirschfelder,  1967)  for  this  classical  assignment  of  electrons  was  introduced  and 
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this  term  will  be  used  in  this  work.  It  is  obvious  that  the  wavefunction  'I'0  is  not  likely 
to  provide  an  adequate  zeroth-order  approximation  for  the  wavefunction,  especially  for 
distances  shorter  than  the  Van  der  Waals  minimum.  In  the  next  sections  two  different 
aspects  of  employing  perturbation  theory  to  intermolecular  forces  will  be  presented.  First, 
I will  briefly  review  the  main  ideas  of  the  Symmetry  Adapted  Perturbation  Theory 
(SAPT). 

Perturbation  Theory  with  Exchange 


In  order  to  present  the  main  versions  of  the  Perturbation  Theory  including  electron 
exchange  the  standard  (or  PA)  theory  may  be  compactly  rewritten  as  follows.  First,  let 
me  introduce  the  reduced  resolvent  operator 


*„  = £ 


Ek  - E° 


(19) 


where  Ek  and  'Pk  are  eigenvalues  and  eigenvectors  of  H°.  Now,  the  Schr0dinger  equation 
reads 


(H°  - E°)$  = (E  - V)tf 


(20) 


Multiplying  this  equation  from  the  left  by  the  resolvent  operator,  using  the  intermediate 
normalization  (^°|^)  = 1 condition  and  noticing  that 


3f?o(tf°-E°)  = 1 - |tf0)($°| 


(21) 


one  finds  easily  that  the  wave  function  can  be  given  as 

= tf°  + »o(E-  V)tf 


(22) 


This  form  suggest  the  iterative  process  to  solve  the  equation.  The  standard  formulae 
of  perturbation  theory  can  be  recovered  by  assuming  that  the  iterative  energy  and 
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wavefunction  converge  and  the  following  expansions  hold  true 

OO 

E = J^E(n) 


n= 1 
oo 

^ ^r(") 
n= 1 


(23) 


Now,  the  m-th  order  energy  and  wavefunction  can  be  calculated  from  recursive  equations 


g(m)  _ 

Tfl 1 /A  J\ 

^r(m)  = EU)^(Tn~j)  ^ 

j= 1 

In  the  PA  method  the  wavefunction  'I'  in  all  orders  of  perturbation  and  in  all  iterations 
is  taken  as  a simple  (or  Hartree)  product  of  two  subsystem  wavefunctions  $a(B)-  Since 
this  approximation  neglects  exchange  a better  choice  of  the  trial  wavefunction  would  be 


(25) 

where  A denotes  the  idempotent  antisymmetrizer  for  all  Na  + Nb  electrons.  The 
antisymmetrizer  can  be  further  decomposed  to 


A = 


Na\Nb\ 
(Na  + Nb)\ 


AaA^(1  + PAB ) 


(26) 


Na(B)  denotes  the  number  of  electrons  of  the  molecule  A(B),  Aa(b>  are  the  antisymmetriz- 

ers  for  A and  B,  while  Pab  is  the  sum  of  permutations  interchanging  at  least  one  pair 

of  electrons  between  A and  B.  This  particular  form  clearly  indicates  the  inter-molecular 

character  of  A.  The  total  Hamiltonian  H commutes  with  A,  but  the  Ho  and  V do  not. 

[H,  A]  = 0 


[tfo,A]^0,[F,A]/0 

A very  interesting  observation  follows  from  these  equations 


(27) 


[#o,  A]  = — [V,  A] 


(28) 
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meaning  that  the  zeroth-order  quantity  (l.h.s.  of  the  last  equation)  equals  the  first-order 
quantity  on  the  r.h.s..  Therefore,  a notion  of  order  is  not  well-defined  in  this  case. 
The  lack  of  well-defined  order  leads  to  fundamental  problems  in  the  case  when  electron 
exchange  cannot  be  neglected.  Unfortunately,  the  function  A$a$b  which  has  proper 
exchange  symmetry  is  not  an  eigenfunction  of  Ho.  There  are  in  principle  two  ways 
of  overcoming  this  difficulty.  One  can  try  to  find  another  zeroth-order  Hamiltonian 
ho  and  perturbation  v such  that  A$a$b  is  an  eigenfunction  of  ho  and  v is  a small 
perturbation.  This  approach  is  called  “symmetric”  perturbation  theory  (Stemheimer,  1954, 
Jansen,  1967).  The  new  zeroth-th  order  Hamiltonian  ho  is  non-hermitian  and  in  practical 
calculations  very  basis  set-dependent.  Despite  many  efforts  the  symmetric  perturbation 
theory  cannot  claim  much  successes  and  therefore  I go  to  the  second  possibility. 

The  second  possibility  consists  of  forcing  the  proper  symmetry  upon  the  wavefunction 
while  keeping  H0  unchanged.  This  Symmetry  Adapted  PT  (SAPT)  has  been  developed 
independently  by  many  authors.  Instead  of  reviewing  them  in  some  order  I will  present 
the  simplest  version  and  an  unified  general  formalism  (Jeziorski  and  Kolos,  1977)  which 
includes  the  majority  of  formulations  as  special  cases. 


Let  us  define  the  properly  antisymmetrized  wavefunction  as  a zeroth-order  trial 
wavefunction 


= iV0  A$J$J  = 


(29) 


Using  this  form  of  the  trial  wavefunction  and  using  the  previously  introduced  iterative 
formalism  one  easily  finds  first-  and  second-order  contributions  to  the  energy 


E(1>  = JV0($°|i/A|$°) 

E<2>  = iV0($° |U»0A(E(1)  - V)|$°) 


(30) 


This  expansion  corresponds  to  the  so  called  MSMA  (Murrell  and  Shaw,  1967,  Musher 


and  Amos,  1967)  formalism.  This  is  the  simplest  form  of  the  SAPT.  Proper  symmetry  is 
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imposed  on  the  zeroth-order  wavefunction  only.  The  simplicity  of  the  MSMA  theory  lies 
in  the  fact  that  no  perturbation  equations  involving  A must  be  solved  in  order  to  calculate 
interaction  energies.  Observing  that  both  A and  the  resolvent  are  hermitian  operators  one 
can  write  the  compact  expression  for  the  second-order  energy 

E(2)  = -V0($^JA(E(1)-V)|$°)  (31) 

where  $pa  denotes  the  wavefunction  obtained  from  the  polarization  approximation,  i.e. 
without  exchange.  It  can  be  proved  (Jeziorski  and  Kolos,  1977)  that  in  all  orders  of 
perturbation  of  the  MSMA  theory  the  antisymmetrizer  A appears  only  in  the  expression 
for  the  energy,  while  the  wavefunction  can  be  calculated  recursively  from  the  equations 
which  do  not  contain  A.  It  can  be  shown  rigorously  that  the  MSMA  methods  provides 
the  correct  asymptotic  expansion  for  the  interaction  energy,  i.e. 

E(n>  - E (p\  = 0(R~k)  (32) 

where  Epa  is  the  energy  obtained  in  the  PA  approximation.  The  difference  between  E 
and  the  PA  part  of  it  in  any  order  of  the  perturbation  theory  may  serve  as  the  definition 
of  the  n-th  order  exchange  energy.  Within  the  MSMA  theory  the  following  splitting  of 
the  n-th  order  energy  holds 

E<”>  = E<”>  + EW  (33) 

The  results  of  the  MSMA  calculations  show  rather  fast  convergence  of  the  method  around 
the  Van  der  Waals  minimum.  It  means  that  for  most  practical  purposes  it  is  possible  to 
approximate  the  total  interaction  energy  as  the  truncated  sum 

AE  a E™  + E<1  + E<?>  + E <1  (34) 

MSMA  can  diverge  at  shorter  distances  and  closer  analysis  (Jeziorski  et  al.,  1980)  reveals 
the  inability  of  the  MSMA  expansion  to  recover  the  total  exchange  contribution  in  a 


finite-order  treatment. 
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As  a remedy  a generalized  iterative  method  was  introduced  (Jeziorski  and  Kolos, 
1977).  The  wavefunction  and  the  energy  are  given 


$ = tf°  + 9?o(E  -V)F* 
E = 


(35) 


The  operators  F and  G are  introduced  to  speed  up  the  convergence  of  the  iterative  process 
and  they  must  satisfy  the  condition 


GV  = F#  = ^ (36) 

where  'P  represents  the  exact  wavefunction.  The  formulae  (35)  provide  an  useful  unifying 
frame  to  analyze  various  realizations  of  the  SAPT.  Assuming  different  forms  of  the 
operators  F and  G and  'I'0  one  can  obtain  different  versions.  For  example,  the  choice  F 
= G = 1 and  = <f>°  corresponds  to  the  standard  Rayleigh-Schrpdinger  perturbation 
theory  or  PA  method.  The  other  choice  is  F = G =1,  but  = jVoA<I>0  This  is  the 
MSMA  method  briefly  described  above.  The  interesting  choice  is  F = A,  G = 1,  and 
'J/°  = 7VoA$o  The  method  is  called  intermediate  symmetry  forcing  (ISF  , Jeziorski  and 
Kolos,  1977,  Jeziorski  et  al.,  1978).  The  convergence  properties  of  the  ISF  expansion 
are  very  good.  The  second-order  energy  displays  the  corect  asymptotic  behavior  at  large 
R identical  to  the  MSMA  method.  However,  the  higher-order  corrections  are  difficult  to 
calculate,  and  they  do  not  have  the  correct  asymptotic  behavior. 

Long-Range  Perturbation  Theory 


At  large  distances  the  exchange  interactions  become  negligibly  small  compared  to  the 
Coulomb  interactions.  Therefore,  the  Hartree  product  becomes  a legitimate 

zeroth-order  wavefunction.  Standard  Rayleigh-Schrpdinger  perturbation  formalism  can 
now  be  applied.  The  formulae  (15),  (16),  (17)  and  (18)  describe  the  first-  and  various 
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second-order  energy  contributions, 
compactly 


The  interaction  operator  V (10)  can  be  written 


v = £ 


sM. 

ra 


(37) 


where  qA  is  the  charge  in  molecule  A,  qB  is  the  charge  in  B and  r denotes  the  distance 
between  the  charges  of  A and  B.  Operator  V can  be  expanded  in  a multipole  series.  This 
is  a very  general  topic  but  here  it  suffices  to  present  the  simplest  approach  based  entirely 
on  classical  electrostatics.  First,  one  has  to  express  r— 1 as  a function  of  the  inverse  of  the 
distance  between  coordinate  origins  for  A and  B.  The  purely  electrostatic  Hamiltonian 
can  be  written  as  the  sum 


v = £?; V = £«?*?  (38) 

» J 

where  <f> i is  the  electrostatic  potential  at  the  charge  i in  molecule  A due  to  the  presence 
of  molecule  B.  The  first  step  in  the  derivation  is  the  Taylor  expansion  of  the  potential 
about  the  coordinate  origin  of  the  molecule  A 

V = T,  « tfV*  + (V«^)or'«  + + ...]  (39) 

i 

The  derivatives  of  the  electrostatic  potential  are  electric  field,  field  gradient  etc.  The 
next  important  simplification  is  the  introduction  of  electric  multipole  tensors.  Many 
general  definitions  are  available  in  the  literature  (Jackson,  1975,  Hirschfelder  et  al.,  1954). 
In  this  work  the  first  few  multipole  moments  are  explicitly  defined  after  Buckingham 
(Buckingham,  1970) 

t*a  = rfr'Q 

i 

®a()  = 2 Xy  — (40) 

I 

An)  (-1)"  A 2n+l  ^ ^ ^ 

n»  Z^«  , drxadrifi"'drxv\vx) 
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Greek  indices  denote  Cartesian  coordinates.  Using  these  definitions  the  expression  for 
V reads 


V = Q-V  - i>.f*  - - 


(41) 


The  Einstein  summation  convention  is  used  in  this  expression.  The  electric  field  at  the 
origin  of  molecule  A due  to  the  charge  distribution  of  B is  denoted  Fa,  while  is 
the  field  gradient.  The  above  defined  electric  moments  are  symmetric  with  respect  to 
interchange  of  any  pair  of  suffixes.  They  also  are  reduced  to  zero  on  contraction  of 
the  tensor  v = 0 . The  second  step  in  the  derivation  of  the  multiple  expansion 

consists  of  expressing  the  potential  and  all  its  derivatives  on  A in  terms  of  the  permanent 
multipoles  on  the  molecule  B.  Again,  the  Taylor  expansion  is  used  and  the  final  result 
for  the  potential  is 

<t>A  = Y,  5*^" 1 “ rjoVaR-1  + l-r)arJ0WaV^R-1  + ...]  (42) 

j 

R denotes  the  distance  between  two  coordinate  origins.  One  can  simplify  the  derivation 
by  introducing  (Jansen,  1957)  the  tensor  T proportional  to  R— <n+1)  and  defined  as  follows 

W.  = (43) 


Since  the  order  of  differentiation  is  irrelevant  the  tensor  T is  symmetric  with  respect  to 
the  exchange  of  any  pair  of  suffices,  it  also  satisfies,  like  the  multipole  moments,  the 
condition  Taapy^_  = 0 . Field  strength  and  field  gradient  can  be  written  in  terms  of 
the  tensor  T as 

R'a  = -QBTq  + iiBTap  - -QByTap1  + —Q’Pft'RaP-rt 
F(*0  = -qBTap  + HyTaPi  - + J^ttB6(Ta0y6e 


(44) 
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To  obtain  the  interaction  Hamiltonian  V one  has  to  combine  the  above  expressions. 
After  some  manipulations  the  final  answer  is 

V = TqAqB  + Ta(qAnB  - qB »£)  + Tap{\qAQBp  + \qBS - pitf) 

36  (45) 

+ Tafa(—qAftap-r  - + 3^®^) 

The  antisymmetry  of  the  last  equation  is  apparent.  The  tensor  T changes  sign  when  the 
vector  R between  two  origin  centers  is  reversed,  or  to  be  more  general 

T'-J  = (46) 

where  subscript  AB  (BA)  explicitly  denotes  the  orientation  of  the  vector  R and  n equals 
the  sum  of  Cartesian  coordinate  indices. 

The  Cartesian  tensors  used  in  the  above  derivation  have  the  advantage  of  rather 
straightforward  algebra.  The  disadvantage  is  that  the  final  expressions  are  highly  redun- 
dant. For  l/2(  n+1  )(  n+2  ) different  components  of  the  symmetrical  tensor  of  rank  n 
only  2n  + 1 contribute  to  the  interaction.  The  geometrical  dependence  of  the  electrostatic 
dependence  is  hidden  in  the  definition  of  the  operator  T.  Again,  it  is  easy  to  present 
a general  expression  but  it  may  be  very  cumbersome  to  work  out  the  particular  case 
of  molecules  of  low  symmetry  in  an  arbitrary  mutual  orientation.  For  these  reasons  an 
alternative  formalism  of  irreducible  tensors  was  introduced.  The  expansion  is  now  given 
in  terms  of  spherical  harmonics  and  Clebsch-Gordan  coefficients.  The  derivation  of  the 
multipole  expansion  (Rose,  1958)  is  however  rather  involved.  The  advantages  of  using 
irreducible  tensors  are  greater  for  higher  order  interactions  and  higher  rank  of  multipoles 
involved.  In  this  work  only  lower  rank  tensors  are  utilized  and  therefore  the  simpler 
Cartesian  tensor  formalism  is  applied. 


CHAPTER  III 


MANY-BODY  THEORY 

In  this  chapter  a brief  summary  of  Coupled-Cluster  (CC)  theory  and  its  finite-  order 
Many-Body  Perturbation  Theory  (MBPT)  approximations  will  be  presented.  It  is  very 
important  to  be  able  to  investigate  the  intermolecular  interactions  with  a theoretical 
method  which  is  size-extensive,  i.e.,  the  energy  of  a system  scales  properly  with  a 
number  of  particles.  The  term  “size-extensivity”  is  used  in  the  sense  introduced  by 
Bartlett  (Bartlett  and  Purvis,  1978,  Bartlett,  1981).  This  term  is  closely  related  to  the 
term  “size-consistency  ” referring  to  a method  whith  a proper  energy  scaling  and  a correct 
reference  function  enabling  a proper  description  of  dissociation  (Pople  et  al.,  1976).  It 
will  be  shown  that  the  CC/MBPT  methods  are  fully  size-extensive. 

The  time-independent  Schrpdinger  equation  for  a molecular  system  is  given  by 


H<H  = E# 


(47) 


where  H is  the  Hamiltonian,  the  wavefunction  is  an  eigenfunction  of  the  operator  H 
and  the  eigenvalue  E is  the  energy  of  the  system.  In  atomic  units,  the  operator  H is 
composed  of  five  parts 


n N V72  n N 7 n N 

"-E^-E^-EE^+E^+E 

.=1  Z A=  1 lMA  »=1  A=1  r'A  i>j  r'>  A>B 


Za%b 

rAB 


(48) 


where  the  capital  labels  A,B,  N refer  to  nuclei  and  lower  case  labels  i,j,n  refer  to  electrons. 
Ma  is  the  mass  of  the  nucleus  A.  The  first  two  terms  are  the  kinetic  energy  terms  for 
electrons  and  nuclei,  respectively.  The  third  sum  is  the  Coulomb  interaction  between 
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electrons  and  nuclei.  The  next  term  represents  the  repulsion  between  electrons,  while  the 
fifth  sum  is  the  repulsion  between  nuclei.  Magnetic  and  relativistic  terms  are  neglected. 
Due  to  the  greater  mass  of  the  nuclei,  they  move  more  slowly  than  electrons.  The 
electrons  may  be  approximated  as  moving  in  a field  of  fixed  nuclei.  This  is  the  essence 
of  the  Bom-Oppenheimer  approximation  which  is  central  to  the  study  of  electronic 
structure.  In  this  approximation  the  kinetic  energy  of  the  nuclei  is  neglected,  while 
the  nuclear  repulsion  term  is  a constant.  The  remaining  terms  of  the  Hamiltonian  define 
the  electronic  problem 


The  electronic  problem  depends  parametrically  on  the  nuclear  coordinates.  The  total 
energy  of  the  system  is  defined  by  adding  the  constant  nuclear  repulsion  term.  The 
Hartree-Fock  method  provides  an  approximate  solution  of  the  electronic  problem. 

Correlated  methods  usually  start  with  some  reference  function  <f>.  In  many  cases 
the  reference  function  is  the  closed-shell  single  determinant  wavefunction  obtained  by 
solving  the  Hartree-Fock  equations.  However,  sometimes  it  is  necessary  to  start  with 
different  reference  function.  The  energy  of  the  reference  function  may  be  expressed  as 


(49) 


The  Schrpdinger  equation  becomes 


He<be  = E e^e 


(50) 


Ere/  = |$>  = (0|ff|0) 


(51) 


The  analysis  of  many-body  methods  can  be  greatly  simplified  if  the  so-called  second 
quantization  or  occupation  number  representation  formalism  is  introduced.  There  exists  a 
one-to-one  correspondence  between  matrix  elements  expressed  in  “ first  quantization”  and 
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second  quantization.  The  Hamiltonian  operator  in  the  occupation  number  representation 
may  be  expressed 

H = Y1  hPiP^  + (wIMpV***  (52) 

p,q 

The  convention  which  will  be  used  in  this  work  is  as  follows.  The  indices 

ij,k,l correspond  to  the  spin  orbitals  occupied  in  the  given  reference  function.  The 

indices  a, b,c,d,.... represent  unoccupied  or  virtual  spin  orbitals.  The  indices  p,q,r,s...are 
unrestricted.  When  the  labels  ij,....a,b,....p,q....are  not  used  as  indices,  they  represent 
creation  (e.g.,  operators  or  annihilation  operators  (e.g.,  a,i,p).  Symbols  hpq  and 

(pq\\rs)  are  defined  as  the  one-and  two-electron  integrals  over  the  basis  functions 

h„  = / - Y, -WW*I  (53) 

J 2 A riA 

and 

(wIM  = J /^(lK(2){^(l-Pi2)}^(l)^(2W^1^2  (54) 

The  variables  Vi  and  o\  are  the  space  and  spin  coordinates  of  electron  number  1, 
respectively.  The  operator  P12  permutes  the  indices  insuring  that  the  integral  <pqllrs> 
is  antisymmetric. 

The  action  of  creation  and  annihilation  operators  can  be  demonstrated  in  the  following 
example 

<^1$)  = | S) 

a'ib'j  |$)  = | D)  (55) 

aUb'jc'k\<f>)  = |T) 

IS>,  ID>,IT>  denote  singly-  , doubly  -,  and  triply-excited  determinant,  respectively.  The 
normal  order  form  is  defined  as  the  product  which  has  all  hole  creation  or  particle 
annihilation  operators  to  the  right  of  the  other  operators.  This  may 

be  accomplished  by  commuting  the  operators.  The  rearrangement  is  done  to  exploit  the 
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fact  that  ^|$)  = 0 a|$)  = 0 i.e.  no  occupied  spin  orbital  can  be  doubly  occupied 

and  an  unoccupied  spin  orbital  cannot  be  annihilated.  In  fact,  these  conditions  illustrate 
the  fundamental  anticommutator  relation 


\p\q]  = pjq  + qp*  = (56) 

If  the  basis  functions  are  non-orthogonal  the  overlap  integral  <plq>  replaces  <$pq. 

The  basic  idea  of  the  cluster  expansion  is  that  the  exact  wavefunction  is  given  by 


l*>  = otp(r)|*) 


(57) 


where  $ is  a known  solution  of  the  independent  particle  problem,  usually  a single  Slater 
determinant.  The  cluster  operator  T is  given  as  the  sum  of  excitations  from  a reference 
function 


where 


T = T1  = T1  + r2  + T>  + ■••• 


Ti  = 

i,a 

T2  = Zt-aUbtj 

i<J 

a<6 

><><* 
a<b  <c 


(58) 


(59) 


The  antisymmetric  t amplitudes  are  to  be  determined.  Inserting  the  exponential  wave- 
function  into  the  Schrpdinger  equation  one  obtains  expression  for  the  energy 


Ecorr  = E - ($\H\$)  = {$\e~THNeT\$)  = ($|(^er)|$)c  (60) 


The  index  C indicates  that  only  connected  diagrams  need  to  be  considered  (Cizek, 
1966,  Bartlett,  1981).  The  suffix  N denotes  that  the  Hamiltonian  is  in  normal  order 
i.e.,  Hn  — H — . The  amplitudes  t are  determined  by  substituting  the 
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CC  wavefunction  (57)  into  Schr0dinger  equation  and  projecting  against  the  excited 
configurations 

($a|(JffiVer)|$)c  = 0 (61) 

In  order  to  evaluate  the  CC  energy  the  cluster  operator  is  expanded.  The  only  nonva- 
nishing terms  are 

Eorr  = {$|  HN(T,  +T2+  Ir,2)|0)c  (62) 

The  truncation  of  the  cluster  operator  in  Eq.  (62)  does  not  depend  on  any  approximate 
version  of  the  CC  theory.  It  rather  reflects  the  fundamental  fact  that  the  Hamiltonian  is 
at  most  a two-electron  operator.  In  actual  calculations  the  correlation  energy  depends  on 
higher  cluster  operators,  i.e.  T3  in  indirect  way,  via  the  amplitude  equations  (61).  In  the 
limit  in  which  all  Tn  cluster  operators  are  include,  the  full  Cl  wavefunction  is  recovered. 
In  practice  the  expansion  of  T is  truncated  leading  to  various  approximations.  Original 
(Cizek,  1969)  formulation  of  the  CC  theory  explicitily  considered  T = T2  or  CCD.  The 
CCSD  model  (Bartlett  and  Purvis,  1978)  is  defined  by  T « T\  + T2  . In  many  cases 
at  least  approximate  consideration  of  triple  excitations  is  of  great  value.  The  family  of 
CCSDT-n  models  (Urban  et  al.,  1985)  can  be  characterized  by  the  approximate  treatment 
of  the  triple  excitations.  The  simplest  model,  CCSDT-1  is  obtained  from  the  following 
truncation  w gTi+T:;  + Th  , subject  to  other  limits  in  the  projection  space. 

A critical  element  in  intermolecular  interaction  theory  is  having  the  correct  results  at 
the  non-interacting  limit.  CC/MBPT  methods  are  extensive,  which  ensure  this  property 
for  interactions  between  closed  shell  molecules.  To  illustrate  this  property,  let  us  consider 
a system  of  localized  noninteracting  electron  pairs  forming  a lattice  of  H2  molecules.  If 
for  simplicity  a set  of  natural  or  Brueckner  orbitals  is  used  to  construct  a reference 
determinant,  only  double-excitations  will  contrbute  to  the  wavefunction.  Obviously,  for 
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N noninteracting  molecules  the  total  energy  is  NEm.  The  exact  wavefunction  for  each 
molecule  i can  be  written 


*(*)  = *o(»)  + x(0  (63) 

where  $o(i)  is  the  Hartree-Fock  determinant  and  x(i)  is  the  sum  of  doubly-excited 
determinants.  The  correlated  wavefunctions  need  not  to  be  orthogonal  (x(*)lx(i))  = 

S6{j  . The  total  wavefunction  for  the  lattice  is 

N 

^ = nW0  + x(*')l  (64) 

t 

It  can  be  readily  noticed  that  the  product  wavefunction  includes  terms  like 

x(0x(i)--x(p)  (65) 

which  correspond  to  simultaneous  double-excitations  on  different  molecules.  However, 
in  the  framework  of  the  Configuration  Interaction  (Cl)  method  these  products  correspond 
to  high-order  simultaneous  excitations. 

The  double-excitation  Cl  wavefunction  can  be  constructed  from  the  reference  function 

*o  = II*o(0 

i 

*DCI  = *0  + *L(k)  (66) 

k^  k 

where  c is  a single  variational  parameter.  The  matrix  elements  are  the  following 

(*I>\hl\xV))  = ? 

(67) 

(<tL(k)\H * |*£(f)>  = + (S  - mhi 

with  a parameter  /?  denoting  correlation  energy  for  an  individual  molecule.  The  DCI 

secular  equation  reads 

AE  — Nc/3  — 0 

(68) 

{/?  - [SAE  + (S'  - l)P}c  = 0 
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where  AE  denotes  the  difference  EL — NEo.  Eq.  (68)  is  a quadratic  equation  for  E.  In 
the  limit  of  large  N the  solution  for  the  energy  is  given  as 


lim  AE  = 

JV— kx> 


(69) 


This  clearly  demonstrates  that  the  DCI  method  cannot  be  size-extensive.  The  correlated 

wavefunction  (64)  can  be  rewritten  as 

N 

41 = nwo + cx(oi = 

N 1=1  <7°) 

n$o(0  + Ec$'  + Ec2*y  + £ «**«»  + •• 

i=l  » i<j  i<j<k 

Eq.  (70)  can  be  further  rearranged  to  an  explicitly  exponential  form.  The  higher- 
excitation  coefficients  are  not  independent  but  are  related  to  the  double-excitations.  In  this 
example  a quadruple  excitation  is  represented  by  the  product  of  two  double  excitations.  In 
the  Coupled  Cluster  theory  the  four-electron  cluster  of  this  type  are  called  disconnected 
clusters.  The  presence  of  the  disconnected  higher-order  clusters  is  essential  for  size- 
extensivity.  It  can  be  also  observed  that  the  standard  Cl  quadruple  excitation  operator 
C4  can  be  expressed  in  terms  of  the  cluster  operators 

C\  = T4  + ir22  + + T1T3  + l-TlT2  (71) 

It  has  been  known  for  a long  time  (Sinanoglu,  1962)  that  the  most  important  part  of  C4 
is  5T2  . The  coupled-cluster  calculations  determine  this  quantity  easily,  while  it  is  very 
difficult  to  introduce  this  part  in  the  Cl  method.  In  other  words,  in  Cl  one  includes  higher- 
excitation  explicitly,  which  extends  the  dimensions  of  the  matrix  problem  enormously, 
while  in  the  Coupled  Cluster  theory  the  most  important  disconnected  parts  of  higher- 
excitations  are  automatically  included.  A more  rigorous  proof  ot  the  size-extensivity  is 
based  on  the  linked-diagram  theorem  (e.g.,  March  et  al.,  1967).  The  main  conclusion  is 
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that  by  eliminating  all  unliked  diagrams  the  resulting  energy  expressions  scale  linearly 
with  the  number  of  particles  and  are  correct  in  the  non-interacting  limit  for  closed  shell 
interactions.  All  finite-order  approximations  to  the  CC  method  known  as  the  Many-Body 
Perturbation  Theory  (MBPT)  (Bartlett,  1981)  are  also  size-extensive  by  virtue  of  the 
same  linked  cluster  theorem. 

It  is  well  known  that  a density  matrix  corresponding  to  an  approximate  wavefunction 
offers  a convenient  way  to  calculate  many  one-electron  properties  (Lpwdin,  1955).  The 
one-electron  reduced  density  matrix  is  defined  as 

7(l,l')=  / *(l,2,...)tf*(l',2 ...)dr2..drN  =£2Wf(1)^(i’)  <72) 

7 P.9 

The  electron  density  is  related  to  the  reduced  density  matrix  by 

MI)  = 7(i,i’)li=i'  = <73> 

P.9 

The  expectation  value  of  any  one-electron  operator  0(1)  is  given  as 

(O)  = J [0(1)M1,  l’)W*7  (74) 

In  the  occupation  number  representation  Dpq  reads 

D„  = <*|o^,|»)  (75) 

When  a Hartree-Fock  wavefunction  for  a closed  shell  system  is  used  as  ^ the  expression 

for  Dpq  becomes  simply  D^F  — Spq  for  p,q  in  occupied  orbitals  and  0 otherwise.  In 
many-body  methods  the  situation  is  more  complicated.  For  example,  the  exponential 
ansatz  leads  to  an  infinite  series  in  Eq.  (72).  In  MBPT  theory  one  can  define  a 
density  by  selecting  a given  MBPT  wavefunction  to  a given  order  in  perturbation.  This 
leads  to  some  problems  even  at  the  first-order.  For  an  SCF  reference  function,  only 
double  excitations  would  contribute  to  the  first-order  wavefunction,  giving  only  double 
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excitations  ('I'(1)|apa9|'I'(1))  contributions  to  the  correlation  correction  for  a property. 
However,  it  is  well  known  that  single  excitations  from  a SCF  reference  function  are 
essential  for  the  correct  description  of  one-electron  properties.  Therefore,  a better 
approximation  to  the  second-order  density  would  be 

D™  = (tfWlaJaJtfW)  + <tf<°>|aj«,l*(2)>  + (¥<2>|ajag|¥(0))  (76) 

where,  besides  the  doubly-excited  function  contributions  from  introduce  single 
excitation  effects  (Amos,  1980).  Similarily,  densities  correct  through  fourth-order  will 
be  considered  in  this  work. 


CHAPTER  IV 


THE  ELECTRIC  POLARIZABILITY  OF  INTERACTING  SYSTEMS 

Introduction 

The  theory  of  intermolecular  interactions  discussed  so  far  deals  with  the  interaction 
energy.  In  this  chapter  changes  in  the  static  dipole  polarizability  caused  by  the  interaction 
will  be  discussed.  Interactions  can  either  modify  the  existing  properties  or  generate  an 
entirely  new  property  like  the  dipole  moment  of  the  He  — Ar  pair.  The  changes  leading 
to  the  occurrence  of  a new  property  are  due  to  the  collision  and  are  often  called  collision- 
induced  phenomena.  The  most  important  are  inductions  of  the  dipole  moment  and  small 
changes  in  the  dipole  polarizability  which  are  responsible  for  the  collision-induced  IR 
and  Raman  specta,  respectively. 

The  translational  absorption  spectrum  of  rare  gas  mixtures  due  to  the  transient  dipole 
moment  induced  in  dissimilar  atoms  by  collisions  was  first  observed  in  the  far  IR  region 
by  Kiss  and  Welsh  (Kiss  and  Welsh,  1959)  and  has  received  a lot  of  attention  since 
then.  The  intensity  of  the  collision-induced  absorption  spectrum  depends  on  the  dipole 
moment  function  /z  (R)  and  the  interatomic  interaction  potential  E (R).  Both  ingredients 
can  be  provided  by  modem  quantum  chemical  computations.  Unlike  the  interatomic 
potential  for  which  many  experimental  sources  of  data  are  available,  \i  (R)  can  be  only 
determined  by  theoretical  calculations.  An  accurate  determination  of  collision-induced 
dipole  moments  encounters  problems  similar  to  the  ones  with  the  Van  der  Waals  potentials 
discussed  before.  The  major  problem  is  the  smallness  of  the  effect.  For  instance,  /z  for 
the  He  — Ar  pair  at  the  collision  diameter  R = 5.65  a.u.  is  estimated  to  be  "0.0075 
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a.u.  ("  0-019  D)  (Meyer  and  Frommhold,  1986).  Perturbational  methods  suffer  from  the 
previously  described  problems  at  short  range,  while  the  SM  calculations  face  the  BSSE 
problems  (e.g.,  Meyer,  1985). 

Intermolecular  interactions  also  give  rise  to  the  Raman  scattering  of  the  homonuclear 
pairs  of  rare  gases.  The  Collision-Induced  Scattering  (CIS)  of  light  effect  was  discovered 
by  McTague  and  Bimbaum  (McTague  and  Bimbaum,  1968).  The  CIS  spectra  are 
generated  by  the  small  incremental  dipole  polarizability  of  the  pair.  The  term  pair 
polarizability  will  be  used  throughout  this  work  for  the  function  a (R) 

a(R)  = aAB(R)  ~ (aA  + QB ) (77) 


where  a denotes  the  dipole  polarizability  of  the  pair  AB  and  individual  atoms.  For  an 
atomic  pair  the  pair  polarizability  is  a diagonal  tensor 


a(R)  = 


<*_L 

0 

0 


where  the  symbols  have  their  usual  meaning, 
trace  a(R)  and  anisotropy  7(R) 


0 

<*_ L 
0 


0 

0 

on 


(78) 


Two  important  quantities  can  be  introduced: 


a(R)  = o(2a-L  + al|) 


(79) 


7 (R)  = <*11  - <*-L 

these  two  quantities  are  experimentally  accesible  in  indirect  ways.  The  trace  function 
is  related  to  the  polarized  collision-induced  Raman  spectra.  The  anisotropy  causes 
depolarization  of  the  spectra  and  the  field-induced  birefrigerence  (Kerr  effect).  For 
densities  comparable  to  the  atmospheric  one  and  for  temperatures  far  from  the  absolute 
zero  the  second  virial  dielectric  coefficient  can  be  related  to  the  trace  according  to  the 
following  formula 


BC(T ) = const. 


+oo 

/ 


a(R) exp( 


0 


(80) 
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E stands  for  the  interaction  energy  and  k denotes  the  Boltzmann  constant.  The  second 
virial  Kerr  coefficient  depends  on  the  anisotropy 

+oo 

Bk(T)  = ^ J -,u(Rha(R)exV(-~)R2dR  (81) 

0 


7a;  refers  to  the  anisotropy  at  a given  optical  frequency  u>,  whereas  70  is  the  static 
anisotropy.  The  expression  for  the  constants  in  (80)  and  (81)  and  more  experimental 
information  can  be  found  e.g.,  in  the  review  by  Frommhold  (Frommhold,  1981).  Most 
experimental  measurments  are  taken  at  frequencies  of  visible  light.  Therefore,  frequency- 
dependent  polarizabilities  should  be  used  to  compare  with  the  experiment.  However,  the 
dynamic  polarizabilities  are  not  readily  available  and  they  differ  from  the  static  values 
(for  the  most  popular  wavelenghts  4880  and  5415  A)  for  the  rare  gase  interactions  by 
less  than  3%  so  the  static  polarizabilities  can  be  rather  safely  used. 

The  simplest,  yet  still  valuable  model  for  the  collision-induced  polarizability  is  the 
classical  Dipole-Induced  Dipole  (DID)  model  (Frommhold,  1981).  In  this  model  an  atom 
(molecule)  is  represented  as  a point  dipole  polarized  by  an  applied  external  electric  field 
Fo  and  the  field  due  to  the  polarization  induced  in  the  partner.  For  a pair  AB  separated 
by  distance  R the  expression  for  the  induced  dipole  moment  can  be  written 

HA  = ao  (-Po  + T{R)pB) 


(82) 


HB  = «0  (-Fo  + T(R)ixa) 


The  suffix  0 denotes  the  isolated  atomic  polarizabilities.  Self-consistent  solutions  for  the 


polarizability  of  A and  B are  given  as 

A qj?(l  +T(R)ag) 

' 1 - ^T(R)aBT(R) 

■ ag(l  + T(R)a#) 

' 1 1 - c,tfT(R)a#T(R) 


(83) 


The  Cartesian  components  of  the  interaction  tensor  operator  T are  defined  in  Chapter 


II.  The  DID  model  holds  only  for  long  range  R and  for  this  reason  it  is  convenient  to 
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express  the  final  forms  for  a and  7 in  terms  of  a series  expansion.  After  assuming  that 


the  atoms  A and  B are  identical  and  truncating  the  series  one  obtains  for  the  trace 


a(R ) « 


1 aj[ 

& 


(84) 


and  for  the  anisotropy 


7(*) 


6a|  6a| 

R*  + R* 


(85) 


At  low  frequencies  the  differences  between  intensities  of  the  depolarized  Raman  spectra 
and  the  DID  model  predictions  are  small.  However,  at  high  frequencies  the  DID 
intensities  can  be  either  3-5  times  too  large  (for  rare  gases,  see  Dacre  and  Frommhold, 
1982)  or  too  small  for  highly  symmetric  systems  CH4,  SF6  (Buckingham  and  Tabisz, 
1978).  The  predictions  for  the  polarized  CIS  spectra  do  not  agree  with  the  experimental 
spectra.  The  DID  models  neglects  long-range  higher  multipoles,  electron  correlation  , as 
well  as  short-range  overlap  and  electron  exchange  effects.  It  should  be  noticed  that  for 
short  distances  R the  formula  (83)  exhibits  a singularity. 

The  first  quantum-mechanical  calculations  of  the  pair  polarizability  were  made  by 
Jansen  and  Mazur  (Jansen  and  Mazur,  1955).  The  corrections  appear  at  the  order  R- 6 
and  are  attributed  to  electron  correlation.  Buckingham  (Buckingham  and  Clarke,  1978) 
developed  a simple  model  that  represents  correlation  effects  as  a change  in  polarizability 
of  one  molecule  in  the  presence  of  another  due  to  the  hyperpolarization  of  the  first  one 
by  the  the  applied  field  and  the  field  gradient  created  by  the  fluctuating  multipoles.  This 
approach  was  extended  (Hunt  et.  al,  1981)  by  introducing  frequency  dependence  of  the 
hyperpolarization  and  polarizability.  These  perturbation  theory-based  formalisms  neglect 
short-range  exchange  and  are  rather  difficult  to  implement.  The  SM  approach  is  not 
free  from  obstacles,  mainly  due  to  basis  set  problems  and  inherent  smallness  of  the  pair 
polarizability.  The  first  SM  calculations  (Buckingham  and  Watts,  1973;  O’Brien  et  al., 
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1973)  for  He2  employed  rather  modest  Gaussian  basis  sets  and  were  performed  at  the 
SCF  level  only.  Using  the  SCF  method  leads  to  purely  repulsive  interatomic  potentials 
for  the  rare  gases  thus  missing  the  important  Van  der  Waals  minimum.  In  the  next  section 
the  calculations  for  Ne2  and  the  new  results  concerning  this  system  will  be  presented. 

Pair  Polarizability  of  Interacting  Ne  Atoms 

Due  to  increasing  experimental  interest  in  collision-induced  phenomena  several 
calculations  of  the  dipole  polarizability  of  the  Ne  — Ne  system  have  been  performed 
(Kress  and  Kozak,  1977;  Dacre  1981;  Dacre  1982;  Dacre  and  Frommhold,  1982).  Kress 
and  Kozak  (KK)  used  an  extended  contracted  Gaussian  basis  set  of  [7s5p3d]  determined 
by  following  the  prescription  by  Werner  and  Meyer  (Werner  and  Meyer,  1976).  The 
calculations  belong  to  the  SM  category.  The  finite  field  Hartree-Fock  (Cohen  and 
Roothaan,  1965)  technique  was  used  with  the  electric  field  strength  of  0.01  a.u.  In 
the  short  range  between  3 and  5 a.u.  these  authors  observed  an  unexpected  shallow 
minimum  on  the  anisotropy  curve.  The  origin  of  this  minimum  was  unclear.  KK 
compared  their  intratomic  potential  to  the  near  HF  quality  results  (Gilbert  and  Wahl, 
1967)  and  found  the  discrepancies  to  be  smaller  than  2%  for  the  very  short  distance  of 
2.5  a.u.  The  observed  minimum  cannot  be  expalained  simply  as  a poor  quality  basis  set 
and  calls  for  another  explanation.  Dacre  (Dacre,  1981)  repeated  KK  calculations.  He 
used  a contracted  [7s4p3d]  basis  set  obtained  from  Huzinaga’s  primitive  (Ils7p/6s3p) 
basis  set  by  augmenting  it  with  ls,lp  and  3d  diffuse  functions.  The  basis  set  is 
similar  to  that  of  KK.  Dacre  computed  the  electron  correlation  by  performing  limited 
CISD  calculations  including  all  single  and  double  excitations  to  a reduced  number  of 
virtual  orbitals.  The  electron  correlation  raised  the  atomic  Ne  polarizability  from  2.3600 
a.u.  to  2.6236  in  good  agreement  with  the  experimental  value  of  2.669a.u.  Since 
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the  limited  Cl  method  is  not  a size-extensive  method  Dacre  had  to  calculate  the  Cl 
polarizability  at  an  interatomic  separation  of  30  a.u.  and  then  subtract  the  average  of 
parallel  and  perpendicular  Cl  polarizabilities.  It  is  obvious  that  such  a procedure  is  only 
an  approximation  and  true  size-extensive  methods  should  be  used  instead.  The  smallness 
of  the  collision-induced  polarizability  can  cause  some  numerical  problems  at  very  long 
distances.  Dacre  investigated  the  possible  effect  of  the  notorious  BSSE  on  the  resulting 
pair  polarizability  at  the  SCF  level  (Dacre,  1981)  and  for  the  CISD  method  (Dacre, 
1982).  It  was  concluded  that  the  BSSE  strongly  affects  the  trace  and  the  valence  only 
Cl  calculations  give  results  very  close  to  the  all-electron  CISD  calculations.  The  short- 
range  anomaly  for  the  anisotropy  could  be  seen  in  Dacre ’s  results,  even  if  he  did  not 
comment  on  this  effect. 

We  use  fully  size-extensive  methods  like  MBPT  and  Coupled  Cluster  theory.  The 
use  of  a method  properly  scaling  with  the  number  of  electrons  is  especially  important 
for  calculating  rather  small  effects  like  the  pair  polarizability,  which  can  be  totally 
hidden  by  empirical  corrections.  The  collision-induced  pair  polarizability  is  obtained 
in  a standard  SM  way:  as  the  difference  of  the  total  polarizability  and  the  sum  of  two 
atomic  polarizabilities.  The  finite  field  method  is  used  to  calculate  the  dipole  electric 
polarizability.  The  strength  of  the  external  electric  field  was  chosen  after  several  test 
calculations.  The  value  0.005  a.u.  does  not  distort  the  electron  density  considerably  and 
still  produces  changes  in  the  total  energy  which  may  be  calculated.  For  a number  of 
selected  points,  a second  field  strength  of  0.010  a.u.  was  used.  The  choice  of  the  basis 
set  is  the  primary  problem  for  calculations  of  second-order  molecular  properties.  Since 
the  basis  used  by  KK  gives  a total  energy  close  to  the  HF  limit  and  at  the  same  time 
yields  an  accurate  value  of  the  polarizability,  a basis  of  similar  quality  was  chosen  for 
our  calculations.  The  basis  consists  of  40  contracted  Gaussian  functions  [7s5p3d]  and 
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was  derived  from  the  uncontracted  (lls6p)  set  (Van  Duijneveldt,  1971).  The  contraction 
coefficients  were  taken  from  benchmark  full  Cl  calculations  (Bauschlicher  et  al.,  1986). 
The  diffuse  orbitals  with  exponents  given  by  Dacre  (Dacre,  1981)  were  added.  All  d 
functions  are  uncontracted.  For  atomic  Ne  calculations  one  f function  with  the  otimized 
exponent  0.28  (Wormer  et  al.,  1983)  was  added.  The  results  did  not  change  much, 
however,  and  the  extra  10  functions  made  the  calculations  excessively  time  consuming. 
Therefore,  most  results  are  obtained  with  a [7s5p3d]  basis  unless  explicitly  denoted 
otherwise.  We  employ  the  approximate  version  of  the  CCSDT  method  in  which  triple 
excitations  are  accounted  for  by  using  converged  values  of  Ti  and  T2  amplitudes  in 
a single  non-iterative  step.  This  method  is  known  as  the  CCSD  +T(CCSD)  (Urban  et 
al.,  1985).  In  contrast  to  CISD  and  other  truncated  Cl-type  expansions  Coupled  Cluster 
theory,  its  approximate  versions  and  MBPT  are  fully  size-extensive,  and,  consequently, 
the  theory  permits  smooth  separation  into  Ne  atoms.  Finite  order  MBPT  results  are 
obtained  in  the  first  few  iterations  of  the  Coupled  Cluster  cycle.  The  results  of  the 
computation  of  a for  Ne  are  presented  in  Table  1.  The  results  obtained  in  these 
calculations  should  be  compared  with  the  experimental  value  of  2.669  a.u.  (Saxon,  1973) 
or  2.671  (in  Taylor  et  al.,  1989).  Our  SCF  results  are  also  close  to  the  near  HF  quality 
results  obtained  by  Taylor  et  al.,  (Taylor  et  al.,  1989)  who  used  a very  extensive  basis 
set  obtained  by  using  an  atomic  natural  orbital  (ANO)  method.  There  is  a discrepancy  of 
0.0088  a.u.  between  our  results  and  the  SCF  results  of  KK.  Adding  one  set  of  f orbitals  to 
the  basis  set  increases  the  results  by  0.001  a.u.  All  results  reported  in  Table  1 are  obtained 
with  two  field  strengths:  0.005  and  0.010  a.u.  The  SCF  energies  were  converged  to  at 
least  10- 11  a.u.  Electron  correlation  amounts  for  13%  of  the  total  dipole  polarizability  as 
given  by  the  full  fourth-order  MBPT  (MBPT-SDQT(4))  computations,  12%  for  the  CCSD 
+T(CCSD)  result  and  10%  at  the  CCSD  level.  The  contribution  of  electron  correlation 
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remains  almost  the  same  when  calculated  with  the  larger  basis  set  [7s5p3dlf]:  13%,  13% 
and  11%  for  MBPT(4),  CCSD  +T(CCSD)  and  CCSD,  respectively. 


Table  1.  Dipole  Polarizability  of  the  Ne  Atom.  All  Values  in  a.u. 


Method 

Basis  [7s5p3d] 

Basis  [7s5p3dlf] 

Ref. 

SCF 

2.361 

2.362 

2.374a,  2.3769b 

MBPT(2) 

2.689 

2.709 

2.7209b 

MBPT(3) 

2.585 

2.600 

2.6065b 

MBPT-SDQ(4) 

2.651 

2.666 

MBPT-SDQT(4) 

2.710 

2.730 

2.7 12a 

CCSD 

2.633 

2.648 

2.61  and  2.64c 

CCSD  -t-T(CCSD) 

2.687 

2.706 

2.70c,  2.71c 

a.  Cemusak,  I.,  Diercksen,  G.H.F.  and  Sadlej,  A.J.,  Phys.  Rev.  A33,  814  (1986) 

b.  McEachran,  R.P.,  Stauffer,  A.D.  and  Greita,  S.,  J.  Phys.  B12 , 3119  (1979) 

c.  Taylor,  P.R.,  Lee,  T.J.,  Rice,  J.E.  and  Almlof,  J.,  Chem.  Phys.  Lett.  163,  359  (1989) 


All  correlated  results  agree  quite  well  with  the  experimental  value.  Only  for  MBPT(3) 
with  the  [7s5p3d]  basis  does  the  discrepancy  exceed  3%.  It  can  be  seen  that  the  full  fourth- 
order  MBPT  calculations  overshoot  the  correct  value.  Enlarging  the  basis  set  by  adding 
an  f function  does  not  change  this  tendency.  The  same  trend  can  be  observed  for  the 
CCSD  + T(CCSD)  results;  better  basis  set  yield  worse  agreement  with  the  experiment. 
Overall  agreement  of  3%  or  better  is  however  acceptable  for  our  purposes. 

In  order  to  assess  the  quality  of  our  basis  set  we  calculated  the  SCF  and  CCSDT-1 
values  of  the  tensors  C,  B and  the  second  hyperpolarizability  7.  The  definitions  of  these 
tensors  will  be  given  in  the  next  chapter.  The  method  of  calculation  will  be  also  discussed 
there.  The  use  of  [7s5p3dlf]  basis  set  yields  the  following  results:  C 1.98  and  2.316, 
B — 12.00  and  — 16.21  and  7 58.5  and  93.4  a.u.  The  first  number  in  each  pair  is  the 
SCF  result,  the  second  is  CCSDT-1/4  value.  These  results  compare  well  with  the  most 
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extended  MBPT(4)  results  (Cemusak  et  al.,  1986):  3.85,  — 17.75  and  104.6  a.u.  for  the 
tensors  C,B  and  7,  respectively. 

Table  2 displays  the  anisotropy  of  Ne2  as  a function  of  R.  The  symbol  D denotes 
double  excitations  only  at  the  second-and  third-order  of  MBPT,  while  SDQ  and  SDQT  are 
shorthand  notations  for  the  MBPT-SDQ(4)  and  the  full  fourth-order  MBPT  calculations, 
respectively.  For  distances  shorter  than  R =4.0  a.u.  Coupled  Cluster  computations  were 
difficult  to  converge  to  the  necessary  limit  (<  10— 7 in  T amplitudes)  hence  only  the  finite 
order  MBPT  results  are  given.  The  results  of  previous  calculations  of  the  anisotropy  7 
(R)  are  given  in  Table  3. 


Table  2 Anisotropy  7(  R ) for  Ne2.  All  Values  in  a.u. 


R 

SCF 

D(2) 

D(3) 

SDQ 

SDQT 

CCSD 

CCSD  + T 

2.5 

0.624 

1.034 

0.869 

1.027 

1.117 

- 

- 

3.0 

0.193 

0.283 

0.249 

0.282 

0.297 

- 

- 

3.5 

0.176 

0.222 

0.213 

0.226 

0.233 

- 

- 

4.0 

0.206 

0.252 

0.243 

0.251 

0.257 

0.253 

0.260 

4.5 

0.216 

0.270 

0.260 

0.269 

0.276 

5.0 

0.204 

0.261 

0.246 

0.256 

0.265 

0.256 

0.267 

6.0 

0.154 

0.201 

0.188 

0.198 

0.207 

0.198 

0.208 

7.0 

0.106 

0.140 

0.129 

0.137 

0.144 

0.136 

0.146 

8.0 

0.072 

0.096 

0.088 

0.094 

0.098 

0.093 

0.098 

10.0 

0.034 

0.046 

0.043 

0.036 

0.046 

0.044 

0.046 

15.0 

0.010 

0.012 

0.012 

0.013 

0.013 

0.012 

0.013 

SCF(68)  denotes  the  use  of  the  basis  [7s4p3d],  while  (80)  denotes  the  basis  set  of  KK. 
The  Cl  calculations  are  of  the  SD  type  with  the  smaller  (68)  basis  set.  The  DF  stands  for 
density  functional  calculations  (Heller  et  al.,  1975).  The  shallow  local  minimum  around 
3.5  a.u.  occurs  in  our  calculations  both  at  the  SCF  and  MBPT  levels.  The  same  behavior 
can  be  seen  in  the  SCF  calculations  by  KK.  Dacre  did  not  report  any  values  for  3.5  and 
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Table  3.  Literature  Values  of  the  Anisotropy  for  Ne2.  All  Data  in  a.u. 


R 

SCF(  68  )a 

CI(68)a 

SCF(  80  )a 

SCF  b 

DF 

2.5 

0.6135 

3.0 

0.1990 

0.2375 

0.1882 

0.1838 

3.5 

0.1670 

4.0 

0.2065 

0.2323 

0.1959 

0.1962 

0.3562 

4.5 

0.2072 

0.3041 

5.0 

0.2054 

0.2340 

0.1961 

0.1962 

0.2480 

6.0 

0.1550 

0.1773 

0.1463 

0.1467 

0.1531 

7.0 

0.1008 

0.0978 

8.0 

0.0732 

0.0833 

0.0691 

0.0693 

0.0656 

10.0 

0.0340 

0.0338 

a)  Dacre,  P.D.,  Can.  J.  Phys.  59,  1439  (1981) 

b)  Kress,  J.W.  and  Kozak,  J.J.,  J.  Chem.  Phys.  66,  4516  (1977) 

c)  Density  functional  calculations,  Heller,  D.F.,  Harris,  R.A.  and  Gelbart,  W.M.,  J. 
Chem.  Phys.  62,  1947  (1975) 

below  3 a.u.,  but  since  his  results  with  basis  (80)  for  3 and  4 a.u.  are  very  close  to  the 
results  of  KK  we  may  conclude  that  his  calculations  also  show  this  unusual  minimum. 
The  difference  of  0.0044  a.u.  between  Dacre’s  and  KK  SCF  result  for  3 a.u.  may  suggest 
that  the  depth  of  the  minimum  is  even  smaller  than  that  reported  by  KK.  From  our  results 
the  depth  of  the  minimum  can  be  estimated  to  be  "0.030  at  the  SCF  level,  the  same  0.030 
for  the  MBPT(2)  and  slightly  less  for  the  full  fourth-order  MBPT:  "0.024  a.u.  The  above 
results  were  obtained  without  any  attempt  to  correct  for  the  basis  set  superposition  error. 
The  discussion  of  this  effect  will  be  given  in  the  next  section.  For  R larger  than  5 a.u. 
the  anisotropy  decreases  monotonically  reaching  the  value  of  "0.010  for  R = 15  a.u.  The 
contribution  of  electron  correlation  to  the  anisotropy  of  Ne2  was  calcualted  previously  be 
Dacre  (Dacre,  1981;  Dacre,  1982)  and  his  results  are  shown  in  Table  3.  For  all  distances 
investigated  we  find  the  correlation  contribution  larger  in  absolute  magnitude  and  more 
significant  in  percentage.  For  R = 3.0  au.  we  calculate  that  the  correlation  amounts  to 
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32%  of  the  total  anisotropy  (0.090  a.u.)  for  MBPT(2),  35%  for  the  full  4-th  order 
(0.104).  The  third-order  contribution  is  the  smallest:  23%.  Dacre  reported  17%  and 
0.0385  a.u.  It  should  be  noticed  that  Dacre’s  Cl  were  obtained  with  a smaller  basis  set 
than  ours.  The  correlation  contribution  is  very  large  for  3 and  3.5  a.u.  then  it  decreases 
to  attain  a minimum  at  4.5  a.u.  and  then  increases  again  to  reach  almost  a constant 
value  from  ~5  to  10  a.u.  The  contribution  becomes  smaller  only  for  large  R = 15.0  a.u. 
The  differences  among  electron  correlation  as  obtained  with  various  many-body  methods 
are  typically  less  than  5%.  For  R = 3.0  the  gap  reaches  12.5%  (the  difference  between 
MBPT(3)  and  the  full  fourth-order).  The  difference  between  second  - and  fourth-order 
is  usually  smaller  than  3%.  CCSD  + T(CCSD)  results  differ  from  the  full  4— th  order  by 
less  than  1%.  Dacre  observed  a similar  behavior  for  the  correlation:  a large  contribution 
at  3 a.u.,  then  a local  minimum  for  R = 4 a.u.,  and  beyond  this  point  rather  a constant 
contribution  up  to  the  distance  of  10  a.u.  What  is  different  however,  is  the  magnitude  of 
the  electron  correlation  he  was  able  to  recover  in  his  CISD  calculations.  His  contribution 
at  the  “constant”  region  is  “13%  in  contrast  to  our  24%  for  MBPT(2)  or  26%  for  the 
CCSD  + T(CCSD)  calculations.  The  density  functional  results  were  obtained  with  an 
assumption  that  the  electric  field-dependent  electron  density  is  nearly  additive,  even  at 
short  distances.  Another  assumption  is  the  use  of  the  energy  functional  of  a homogeneous 
electron  gas  at  each  point  in  space.  The  results  are  much  too  high  at  short  distances  and 
begin  to  agree  with  the  MO  calculations  only  for  distances  longer  than  6 a.u. 

Table  4 dispalys  the  trace  a (R)  as  obtained  by  different  MBPT  and  Coupled  Cluster 
methods.  For  comparison  previous  results  are  shown  in  Table  5.  Like  for  the  anisotropy, 
the  SCF  and  Cl  results  of  Dacre,  SCF  results  of  KK  and  density  functional  computation 
results  are  presented.  The  same  shorthand  notation  as  in  Table  2 is  used  here.  The  SCF 
values  agree  well  with  the  results  of  Dacre,  the  difference  varying  from  0.005  a.u.  at 
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R equal  to  3.0  a.u.,  0.004  at  5.0  a.u  and  0.0024  for  R = 6.00  a.u.  For  R longer  than 
6 a.u.  the  trace  becomes  very  small  and  the  numerical  errors  make  the  results  oscillate 
around  zero.  Like  in  the  case  of  the  anisotropy  the  correlation  contribution  calculated 
with  MBPT/Coupled  Cluster  methods  is  almost  two  times  larger  than  Dacre’s  CISD 
results.  For  R = 2.5  a.u.  the  correlation  contribution  makes  the  trace  less  negative,  but 
it  can  be  seen  that  the  differences  among  various  correlated  results  are  large:  MBPT(3) 
gives  — 0.213  while  MBPT-SDQT(4)  yields  only  — 0.133  a.u.  For  all  other  short  range 
distances  the  differences  among  correlated  results  are  smaller  and  show  the  expected 
pattern.  The  correlation  contribution,  as  given  by  fourth-order  MBPT,  increases  from 
13%  at  3.0  a.u.  to  21%  for  R equal  5 a.u.  Due  to  the  above  mentioned  convergence 
problems  the  Coupled  Cluster  results  are  not  available  for  all  values  of  the  interatomic 
distance.  When  available  the  CCSD  results  differ  from  the  MBPT(4):  e.g.  at  R = 5.0  a.u. 
the  CCSD  + T(CCSD)  calculations  yield  — 0.014  a.u.  which  is  half  the  MBPT  result. 
It  is  not  clear  whether  these  differences  reflect  the  importance  of  higher  excitations  or 
rather  numerical  errors.  For  R greater  than  6.0  a.u.  the  Coupled  Cluster  results  and  the 
MBPT  oscillate. 

From  Table  5,  it  can  be  seen  that  the  similar  oscillatory  results  occur  for  the  CISD 
calculations.  Our  SCF  results  are  very  close  to  the  results  of  KK  and  only  for  distances 
longer  than  6.0  a.u.  do  differences  emerge.  For  distances  longer  than  4.0  a.u.  the 
density  functional  calculations  yield  more  acceptable  results  than  for  the  anisotropy. 
The  asymptotic  behavior  of  the  density  functional  trace  is  wrong,  but  all  other  SM- 
type  calculations  suffer  from  numerical  problems  in  this  region.  Inclusion  of  electron 
correlation  does  not  make  the  results  more  numerically  stable.  The  use  of  the  non-size 
extensive  CISD  method  is  even  more  problematic. 
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Table  4.  Trace  a (R)  for  Ne2-  All  Values  in  a.u. 


R 

SCF 

D(2) 

D(3) 

SDQ 

SDQT 

CCSD 

CCSD  + 
T 

2.5 

-0.258 

-0.191 

-0.213 

-0.157 

-0.133 

3.0 

-0.289 

-0.336 

-0.315 

-0.323 

-0.332 

3.5 

-0.193 

-0.237 

-0.215 

-0.226 

-0.235 

4.0 

-0.108 

-0.136 

-0.123 

-0.131 

-0.137 

-0.127 

-0.132 

4.5 

-0.053 

-0.066 

-0.060 

-0.064 

-0.067 

5.0 

-0.022 

-0.027 

-0.028 

-0.028 

-0.028 

-0.020 

-0.014 

6.0 

-0.002 

0.001 

-0.008 

0.004 

0.005 

0.000 

0.001 

7.0 

0.003 

0.004 

-0.003 

0.002 

0.002 

0.001 

0.002 

8.0 

0.002 

0.004 

-0.000 

0.003 

0.005 

0.004 

0.005 

Table  5.  Trace  of  Ne2.  Literature  Values  in  a.u. 


R 

SCF(68)a 

SCF(80)a 

CI(68)a 

SCF  b 

DF 

2.5 

-0.2751 

3.0 

-0.2874 

-0.2945 

-0.3157 

-0.2949 

3.5 

-0.1960 

4.0 

-0.1061 

-0.1104 

-0.1232 

-0.1104 

-0.2023 

4.5 

-0.0553 

-0.0802 

5.0 

-0.0207 

-0.0247 

-0.0271 

-0.0247 

-0.0286 

6.0 

0.0002 

-0.0026 

-0.0028 

-0.0027 

-0.0010 

7.0 

0.0010 

0.0010 

8.0 

0.0027 

0.0013 

0.0009 

0.00013 

0.0004 

a)  Dacre,  P.D.,  Can.  J.Phys.  59,  1439  (1981) 

b)  Kress,  J.W.  and  Kozak,  J.J.,  J.  Chem.  Phys.  66,  4516  (1977) 

c)  Heller,  D.H.,  Harris,  R.A.  and  Gelbart,  W.M.,  J.  Chem.  Phys.  62,  1947  (1975) 


Secondary  Basis  Set  Superposition  Error 


All  SM  results  are  contaminated  by  the  systematic  Basis  Set  Superposition  Error. 
BSSE  affects  electric  properties  of  the  interacting  molecules  as  well.  In  Chapter  II 
the  most  important  methods  designed  to  eliminate  or  reduce  this  error  were  presented. 
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The  general  conclusion  is  that  despite  much  effort  the  elimination  of  the  BSSE  is  still 
a debatable  matter.  For  molecular  properties  other  than  energy  the  basis  superposition 
problem  is  even  less  known.  Karlstrpm  and  Sadlej  (Karlstrpm  and  Sadlej,  1982)  analyzed 
this  problem  for  the  (H 20)2  system.  The  electron  density  and  the  polarizability  may 
change  considerably  from  the  isolated  molecules  to  the  SM  (dimer)  basis  set.  This  is 
especially  true  for  the  minimum  and  small  basis  sets.  These  changes  are  not  accounted  for 
by  the  standard  Boys-Bemardi  or  CP  correction  and  can  be  called  higher-order  basis  set 
superposition  errors.  The  secondary  BSSE  can  be  formally  defined  for  a given  electric 
property,  e.g.,  dipole  moment  as 

&HA  = U B)  - n(A)  (86) 

where  the  notation  concerning  the  basis  sets  employed  is  the  same  as  in  Chapter  II.  How- 
ever, contrary  to  the  BSSE  for  energy  the  higher-order  or  secondary  superposition  effects 
may  lead  to  improvment  of  the  calculated  interaction  energies  through  the  improvement 
of  the  electric  moments  of  the  molecules.  It  should  be  noted  that  the  moments  can  be 
either  increased  or  decreased  by  the  basis  set  superposition.  One  can  anticipate  that  the 
electrostatic  part  of  the  interaction  energy  can  be  better  described  when  the  superposition 
of  the  basis  is  present.  Superposition  of  the  basis  leads  also  to  changes  in  the  polarizabil- 
ities. Larger  basis  sets  almost  always  lead  to  an  increase  in  the  polarizabilities  and  thus 
change,  considerably,  the  induction  part  of  the  interaction  energy.  For  the  discussion  of 
the  secondary  BSSE  on  the  interaction  energies  Karlstrpm  and  Sadlej  introduced  ratios  of 
the  dipole-dipole  and  induction  energies  calculated  with  different  sets  of  water  molecule 
dipole  moments  and  polarizabilities.  The  minimal  dimer  basis  set,  minimal  basis  set  for 
isolated  molecule  and  near-HF  values  were  analyzed  in  this  study.  The  main  conclusion 
was  that  both  electrostatic  and  induction  interaction  energies  enhance  considerably  in  the 
regions  close  to  the  SCF  energy  minima.  The  change  of  the  induction  energy  is  much 
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faster  than  the  electrostatic  energy.  A minimal  basis  set  overestimates  the  electrostatic 
part  while  seriously  underestimating  the  induction  energy.  The  compensation  of  these 
two  effects  does  not  occur  for  the  water  dimer  and  cannot  be  expected  to  hold  in  general. 
Scheiner  and  his  group  (Latajka  and  Scheiner,  1987a;  Latajka  and  Scheiner,  1987b)  ana- 
lyzed the  secondary  BSSE  on  the  SCF  and  MBPT(2)  level  of  theory.  He  recommended 
the  subtraction  of  both  primary  and  secondary  BSSE.  He  also  found  that  the  secondary 
superposition  effect  is  highly  anisotropic  and  this  can  affect  the  interaction  energy.  Sec- 
ondary BSSE  was  studied  in  the  IR  and  Raman  intensity  for  the  (HF)2  dimer  (Szczesniak 
and  Scheiner,  1986).  It  was  concluded  that  the  IR  intensities  are  relatively  superposition 
error-free  but  the  Raman  intensity  is  rather  sensitive  to  this  error.  Due  to  the  large  error 
the  use  of  the  popular  basis  set  of  the  6-3 1G  etc.  family  for  Raman  intensity  studies 
is  not  recommended.  The  latter  conclusion  was  made  on  the  basis  of  SCF  calculations 
employing  small  to  modest  basis  functions.  It  should  also  be  noticed  that  the  Raman 
intensity  depends  on  the  derivative  of  the  polarizability,  so  the  large  superposition  error 
may  not  entirely  come  from  the  polarizability  itself. 

We  accounted  for  the  basis  set  superposition  error  by  using  the  standard  counterpoise 
correction  (CP).  The  CP  method  results  in  obtaining  parallel  and  perpendicular  com- 
ponents of  the  atomic  dipole  polarizability  as  a function  of  the  interatomic  distance  R. 
In  Table  6,  we  present  the  perpendicular  component  of  the  polarizability  for  a selected 
number  of  points.  Differences  at  the  SCF  level  are  very  small:  0.0014  a.u.  from  3.0 
to  7.0  a.u.  Differences  at  the  correlated  level  are  slightly  larger.  MBPT(2)  yields  the 
difference  of  0.0029  (for  R = 4.0  a.u.),  while  the  full  fourth-order  gives  0.0032  a.u.  the 
SCF  result  reaches  a constant  value  of  2.3605  at  R = 6.0  a.u.,  while  the  perpendicular 
polarizability  at  the  MBPT  level  continues  to  decrease  slightly.  Concluding  we  may  say 
that  the  BSSE  effects  for  the  perpendicular  component  are  below  the  acceptable  accuracy 
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of  our  computations  and  hence  can  be  neglected,  at  least  for  distances  longer  than  4.0 
a.u.  Table  7 displays  the  CP  corrected  values  of  the  parallel  component  of  the  atomic 


polarizability  tensor. 

Table  6.  Counterpoised  Corected  Ne  Perpendicular  Polarizability  in  a.u. 


R 

SCF 

D(2) 

D(3) 

SDQ 

SDQT 

3.0 

2.3619 

2.6945 

2.5895 

2.6563 

2.7167 

3.5 

2.3618 

2.6936 

2.5892 

2.6559 

2.7161 

4.0 

2.3616 

2.6920 

2.5875 

2.6539 

2.7138 

4.5 

2.3616 

2.6920 

2.5872 

2.6544 

2.7136 

5.0 

2.3608 

2.6904 

2.5864 

2.6520 

2.7120 

6.0 

2.3605 

2.6895 

2.5860 

2.6516 

2.7108 

7.0 

2.3605 

2.6891 

2.5858 

2.6514 

2.7106 

Table  7.  Counterpoised  Corected  Ne  Parallel  Polarizability  in  a.u. 


R 

SCF 

D(2) 

D(3) 

SDQ 

SDQT 

3.0 

2.3657 

2.7070 

2.5986 

2.6681 

2.7308 

3.5 

2.3654 

2.7036 

2.5963 

2.6653 

2.7273 

4.0 

2.3653 

2.7015 

2.5950 

2.6635 

2.7250 

4.5 

2.3652 

2.7012 

2.5944 

2.6632 

2.7240 

5.0 

2.3656 

2.7001 

2.5943 

2.6620 

2.7229 

6.0 

2.3658 

2.6989 

2.5929 

2.6602 

2.7207 

7.0 

2.3649 

2.6969 

2.5920 

2.6589 

2.7190 

The  influence  of  the  basis  set  superposition  effect  on  the  parallel  component  of  the 
polarizability  is  stronger.  Generally,  the  values  of  the  polarizability  are  larger  now. 
At  the  interatomic  distance  equal  to  4.0  a.u.  the  difference  between  CP-corrected  and 
a “raw”  polarizability  as  calculated  by  the  SCF  method  is  0.0047  a.u.  about  6 times 
greater  than  for  the  perpendicular  component.  For  MBPT(2)  the  change  is  0.0120  and 
0.0148  a.u.  for  the  MBPT-SDQT(4)  calculations.  Like  for  the  perpendicular  component 
the  correlated  results  decrease  smoothly  with  increasing  distance,  only  at  the  SCF  level 


48 


can  very  small  irregularities  be  seen  for  R equal  to  5.0  and  6.0  a.u.  From  the  above  data 
the  CP-corrected  anisotropy  and  trace  can  be  calculated.  The  results  for  the  CP-corrected 
anisotropy  function  7CP  (R)  are  given  in  Table  8.  Table  9 contains  the  CP-corrected 


values  of  the  tarce  a (R). 

Table  8.  Counterpoised  Corrected  Anisotropy  of  Ne2  in  a.u. 


R 

SCF 

D(2) 

D(3) 

SDQ 

SDQT 

SCFa 

3.0 

0.186 

0.258 

0.230 

0.259 

0.269 

0.1845 

3.5 

0.169 

0.202 

0.199 

0.208 

0.210 

4.0 

0.198 

0.233 

0.229 

0.233 

0.234 

0.1987 

4.5 

0.216 

0.252 

0.246 

0.251 

0.256 

5.0 

0.196 

0.241 

0.231 

0.237 

0.243 

0.1986 

6.0 

0.143 

0.182 

0.174 

0.181 

0.187 

0.1438 

a)  SCF(80)  calculations,  Dacre,  P.D.,  Can.  J.  Phys.  59,  1439  (1981) 


Table  9.  Counterpoised  Corrected  Trace  of  Ne2  in  a.u. 


R 

SCF 

D(2) 

D(3) 

SDQ 

SDQT 

SCFa 

3.0 

-0.294 

-0.356 

-0.328 

-0.341 

-0.354 

-0.2965 

3.5 

-0.197 

-0.253 

-0.226 

-0.242 

-0.254 

4.0 

-0.112 

-0.149 

-0.132 

-0.142 

-0.151 

-0.1121 

4.5 

-0.057 

-0.078 

-0.068 

-0.076 

-0.078 

5.0 

-0.027 

-0.036 

-0.036 

-0.037 

-0.039 

-0.0266 

6.0 

-0.005 

-0.007 

-0.013 

-0.003 

-0.003 

-0.0042 

a)  See  Table  8 


As  can  be  seen  from  Table  8,  the  anisotropy  with  the  counterpoise  corrections  behaves 
similarily  to  the  “raw”  anisotropy.  At  distances  between  3 and  5 a.u.  the  local  minimum 
occurs  and  its  depth  is  estimated  to  be  ~0.05  a.u.  both  at  the  SCF  and  correlated  levels 
of  computations.  The  effect  of  the  CP  correction  on  the  anisotropy  is  to  reduce  the 
calculated  value.  The  CP  corrections  applied  at  the  MBPT  level  of  theory  are  larger 
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than  the  corresponding  SCF  CP  corrections.  For  example,  at  R = 4.0  a.u.  the  SCF  CP 
correction  is  0.008  a.u.  while  the  second-order  CP  correction  is  0.019  and  0.023  a.u. 
for  the  full  fourth-order.  Dacre  (Dacre,  1981)  also  investigated  the  basis  superposition 
effect  on  the  SCF  results.  The  results  obtained  with  a larger  basis  containing  80  functions 
are  given  for  comparison  in  the  last  column  of  Table  8.  In  general  both  SCF(68)  and 
SCF(80)  results  after  corrections  agree  better  among  themselves  and  are  in  very  good 
agreement  with  our  uncorrelated  results.  The  basis  extension  corrections  reduce  the  trace 
value  and  the  trace  becomes  more  negative.  As  for  the  uncorrected  case  for  R greater 
than  6 the  results  become  less  reliable  due  to  numerical  instabilities.  For  R equal  6 a.u. 
the  counterpoise  corrected  trace  are  all  negative  but  the  differences  among  correlated 
methods  are  still  large. 

It  may  be  interesting  to  compare  the  secondary  basis  set  superposition  to  the  standard 
BSSE.  Table  10  shows  the  interaction  energy  calculated  with  and  without  the  Boys- 
Bemardi  correction  for  three  intemuclear  distances. 


Table  10.  Interaction  Energy  of  Ne2.  Energies  in  //Hartree,  Distance  in  a.u. 


R 

SCF 

D(2) 

SDQT 

CCSD 

CCSD  + T 

5.0 

713.7 

303.6 

218.4 

243.6 

195.5 

5.0  CP 

739.1 

516.0 

448.4 

457.2 

421.5 

6.0 

43.9 

-171.6 

-206.8 

-187.7 

-211.7 

6.0  CP 

63.1 

-48.4 

-75.8 

-64.7 

-82.6 

7.0 

-9.3 

-119.5 

-134.4 

-124.1 

-135.4 

7.0  CP 

5.3 

-43.4 

-54.8 

-48.8 

-56.8 

The  BSSE  change  the  interaction  drastically.  Notice  that  without  the  CP  correction 
there  is  an  artificial  minimum  at  the  SCF  level  for  R equal  to  7.0  a.u.  Since,  the  SCF 
method  does  not  include  dispersion  interactions  responsible  for  any  long-range  minima 
the  minimum  of  — 9.3  ^Hartree  must  be  due  to  the  basis  set  expansion  error.  The 
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correlated  results  are  even  more  susceptibile  to  the  BSSE.  For  R = 5.0  a.u.  the  BSSE 
inreases  the  interaction  energy  making  the  potential  curve  more  repulsive.  The  change 
is  of  the  magnitude  of  the  repulsion  itself.  For  the  attractive  region  the  CP-corrected 
energies  are  smaller,  hence  the  potential  well  is  less  attractive.  It  is  obvious  that  use  of 
the  uncorrected  potential  gives  a qualitatively  wrong,  unphysical  answer.  On  the  other 
hand,  the  effect  of  the  basis  extension  upon  the  dipole  polarizability  is  shown  to  be  small 
for  the  anisotropy  and  modest  for  the  trace.  In  neither  case  does  the  secondary  BSSE 
change  the  results  substantially. 


CHAPTER  V 


ELECTRIC  MULTIPOLES,  POLARIZABILITY  AND  HYPERPOLARIZABILITY 

Charge  Perturbation  Method 


Electrostatic  or  Coulombic  and  induction  contributions  to  the  interaction  energy  can 
be  evaluated  using  classical  electrostatic  concepts.  For  a great  variety  of  systems  these 
interactions  play  a dominant  role.  The  Coulombic  interactions  can  be  represented  as  the 
interactions  between  two  sets  of  permanent  electric  multipoles  located  on  molecules  A 
and  B.  The  multipoles  on  A create  an  inhomogenous  electric  field,  which  distorts  the 
charge  distribution  of  the  molecule  B giving  rise  to  the  induction  energy.  The  ability  of 
molecule  B to  be  distorted  by  the  field  caused  by  A is  measured  by  the  polarizability 
tensors  of  B (Buckingham,  1967).  Therefore,  in  the  long-range  limit  the  induction  can 
be  expressed  as  a sum  involving  dipole,  quadrupole  and  other  static  polarizabilities  of 
B subject  to  the  Taylor  series  expansion.  The  expansion  introduces  the  gradients  and 
higher  derivatives  of  the  field  at  A.  Coulombic  and  induction  interactions  can  then  be 
expressed  as 


Ecu,  = T,V  + T „(«X  - 


eLi  = - ^F?F?F?  - ^J^F^  - ^C^F&FB 

The  tensor  T was  already  introduced  in  Chapter  II.  Other  symbols  will  be  explained 
below.  It  is  obvious  that  the  problem  of  calculating  the  Coulombic  and  induction  energy 
revolves  around  the  accurate  and  efficient  calculation  of  the  permanent  moments  and 
polarizabilities.  The  multipole  moments  are  first-order  properties  and  can  be  calculated 


(87) 
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as  the  first  derivatives  in  the  Taylor  series  expansion  of  the  energy  or,  from  an  expectation 
value-type  expression.  The  former  was  employed  in  the  Ne2  study,  while  the  latter  will  be 
considered  in  the  following  section.  The  equivalence  between  the  two  expansions  implies 
the  satisfaction  of  the  Hellmann-Feynman  theorem.  The  relationship  between  these  two 
different  approaches  will  be  discussed  in  the  next  section.  The  polarizabilities  are  second- 
or  higher-order  properties  and  are  more  difficult  to  calculate.  One  can  use  either  analytical 
differentation  methods  (e.g.  Amos,  1987)  or  employ  more  straightforward  finite-field 
techniques.  In  this  work  the  modification  of  the  finite-field  technique  called  the  charge 
perturbation  method  (McLean  and  Yoshimine,  1967a,  1967b)  is  used. 

In  the  finite-field  technique  a molecule  (or  atom)  is  placed  in  a weak  external 
homogeneous  electric  field  and  then  the  total  molecular  energy  is  calculated.  In  principle 
the  first  -,  second  -,  and  higher-order  electric  properties  of  the  system  can  be  calculated 
by  numerical  differentation  of  the  total  perturbed  energy.  The  pertinent  formulae  can 
be  found  in  the  literature  (Bartlett  and  Purvis,  1979).  Since  the  only  property  needed 
to  determine  the  polarizability  in  this  method  is  the  total  pertubed  energy  there  is  no 
fundamental  difficulty  in  using  correlated  methods  like  MBPT/CC.  In  this  chapter  all 
results  obtained  by  numerical  differentiation  will  be  denoted  as  finite-field  (FF)  results, 
however,  it  is  necessary  to  distinguish  between  derivatives  of  the  energy  and  of  a moment. 
The  normal  finite-field  method  simply  adds  to  the  one-particle  Hamiltonian  an  electric 
field  operator  multiplied  by  a small  constant,  A;  • r;  where  A represents  the  field 

t 

strength  in  the  appropriate  direction.  An  alternative  to  this  procedure  is  to  generate  the 
electric  field  by  an  extra  charge  placed  at  a distance  R from  the  molecule.  The  field  is  not 
homogeneous  and  the  electrostatic  potential  can  be  expanded  in  a Taylor  series  around 
the  chosen  point  0 inside  the  molecular  charge  distribution 

cf>(x,  y,  z)  = <j>{ 0)  + ra^( 0)  + ^rar^"afi(0)  + ^rarjjr7^7(0)  + ...  (88) 
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Primes  denote  partial  derivatives  of  the  potential.  The  derivatives  are  totally  symmetric 
and  are  related  to  the  electric  field  strength  and  field  gradient  through  the  following 
expressions 

(89) 


'—(£). 


Throughout  this  work  the  derivatives  of  the  electrostatic  potential  are  always  calculated 
at  the  center  of  the  mass.  The  total  Hamiltonian  of  the  system  can  be  written  as  the  sum 


flel  At  Nnucl  ry 

hw.*)=h<“’-E^+E  ZlQ 


(90) 


j- i lri  - ^1  jTj  I rJ  ~ R\ 

Expanding  the  Hamiltonian  about  the  unperturbed  H(0)  in  terms  of  the  molecular  multipole 
moments  one  obtains  the  following  expression: 


H = H(0)  - HaFa  - - 9apFa/ $ - — tta/3-fFapy  - ... 


(91) 


while  the  corresponding  total  perturbed  energy  E reads 

E = E(0)  - HaFa  - -QapFap  - — SlaPyFapy  - ... 

2 ®apFaFp  -Aatp^FaFpy  ~ Capty^FapFyg  ..  (92) 

g/^a/?7  FaFpF^  — —Ba  p^fiFaFpFyg  — ■^^/yapy^FaFpFyF^  .. 

//,  0,  ft  denote  the  dipole,  quadrupole  and  octopole  moments  of  the  unperturbed  molecule. 
a,  /3,  7 are  the  dipole  polarizability,  hyperpolarizability  and  second  dipole  hyperpolariz- 
ability respectively.  Dipole-quadrupole,  quadrupole-quadrupole,  dipole-octopole  tensors 
are  denoted  by  the  symbols  A,  C and  E,  respectively.  The  third-order  tensor  B is  the  mixed 
dipole-dipole-quadrupole  hyperpolarizabilty.  In  Eq.  (92)  the  first  line  corrsponds  to  the 
first-order,  the  second  to  the  second-order  energy  and  the  last  line  contains  the  third  - and 
fourth-order  (7)  terms.  The  same  Taylor  expansion  can  be  applied  to  the  external  field- 
dependent  multipole  moments  of  the  perturbed  molecule.  The  coefficients  in  the  power 
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series  are  the  permanent  (field-independent)  mutipole  moments  and  polarizabilities 

Hai.F yi  F ar/?)  ••)  = T aa/3F^  "h  ^Aa^yF^y  "I-  ~^Fa^y^F^Fy^ 

Eap,  ••)  ==  ©a/J  “h  Ay^pFy  *H  Ca/}tf$Fyg  + —ByfiapFyFfj  (93) 

^aP-f{For)  Fafii  ••)  = ^a/?7  T"  ^6,afiyF § + ~H^ea^yF^c 

The  new  symbol  H corresponds  to  the  quadrupole-octopole  polarizability,  a second-order 
tensor.  The  perturbed  or  field-dependent  moments  must  be  calculated,  while  the  field 
parameters  can  be  easily  obtained  for  a given  charge  Q and  a distance  R.  In  practice  one 
can  use  symmetry  arguments  to  simplify  the  geometrical  dependence  of  the  moments. 
Tensors  aap,  pa ^ , 7a/?76  and  moments  Qap,  are  unchanged  on  permuting 

the  indices.  The  tensor  Aa  p^  is  symmetric  in  the  indices  /3j,  Cap  ^ is  symmetric 
in  either  pair  of  indices,  while  it  holds  that  Captl6  = Cy6<ap  . The  tensor  Ea  p $ is 
symmetric  under  permutations  of  the  indices  (3^6 . The  tensor  Bap  ^ does  not  change 
under  permutations  within  either  pair,  but  unlike  C,  the  exchange  property  does  not  hold. 
Other  useful  relationships  are  a consequence  of  the  fact  that  Laplace’s  equation  is  satisfied 
inside  the  charge  distribution  (molecule).  It  can  be  shown  (Buckingham,  1978)  that  the 
following  relationships  are  satisfied:  Oaa  = Qap p = Aa  pp  = 0. 

The  unknown  tensors  a,  /?,  A... .can  be  obtained  by  solving  the  set  of  equations  (92) 
and  (93).  For  a spherical  system  like  the  fluoride  anion  F- (*S)  in  the  presence 
of  a single  external  charge  Q at  a distance  R from  the  origin,  the  problem  of  the 
polarizability  tensor  determination  can  be  readily  solved.  First,  all  multipole  moments 
Va  = = 0 and  (3apy  = Aajy  = Ea^ = 0 . The  reamining  tensors  are 

isotropic.  The  electric  field  created  by  the  charge  Q evaluated  at  the  nucleus  is  — Q/R2. 
Equations  (92)  then  reads 


E (Q,  R ) = E<°)  - ^ 


3 CQ2  BQ 3 
2 i?6  + 2 R7 


7 Q4 
24  i?8 


(94) 
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Equations  for  the  field-induced  dipole  and  quadrupole  are  similarly  very  simple 

m <*Q  , BQ2  703 

M,R)~  r2+  rS  r6 

am  m 3C<?  i bq2 
0(Q,R)-  r3  + Rt 

The  equation  for  the  polarizabilities  are  now  given  as: 

R 2 


(95) 


a = 


12  Q 
R 6 

7 “ 2Q* 

R 3 


M-2Q)  - K2Q)  ~ M-Q)  + 8 M)} 

{fi(-2Q)  - ti{2 Q)  - 2 fi(-Q)  + 2n(Q)} 


(96) 


C = ~6Q  ‘ {Q{~Q)  ~ Q(Q)} 

B=~-  {©(-<?)  + 0(Q)} 

As  in  the  standard  FF  method  the  charge  perturbation  requires  several  calculations 
to  determine  the  desired  polarizability  tensor.  The  external  charge  Q has  a similar 
significance  to  the  electric  field  strength  in  the  standard  FF.  There  is  no  method  to 
choose  the  value  of  Q (and  R)  but  trial  and  error.  In  order  to  obtain  reliable  results, 
extrapolation  to  zero  must  be  performed. 

Before  proceeding  further  some  comments  on  the  above  derivation  are  in  order. 
The  multipole  moments  used  above  are  traceless  tensors.  In  other  words  the  general 
definition  of  the  mutipole  moment 


«<n)  = i J p(r)rndv  (97) 

V 

where  p is  the  electric  charge  density,  r"  denotes  the  direct  product  of  the  n factors  r and 
the  integration  is  carried  out  over  the  whole  volume  containing  the  charge  density,  has  to 
be  modified.  The  symmetry  relations  following  from  Laplace’s  equation  are  used  and  the 
detailed  derivation  can  be  found  elsewhere  (McLean  and  Yoshimine,  1967b).  For  higher- 
order  moments  the  Cartesian  tensors  soon  become  redundant  and  rather  cumbersome  to 
use.  The  use  of  the  spherical  tensor  formulation  for  multipole  moments  (Gray  and  Lo, 
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1976)  is  then  recommended.  A new  reformulation  in  terms  of  the  unabridged  Cartesian 
tensors  has  been  presented  recently  (Applequist,  1984,  1985).  The  new  formalism  leads 
to  simpler  relationships  between  derivatives  of  the  energy  with  respect  to  the  field  and 
the  multipole  moments  and  polarizability. 

It  should  be  recognized  that  the  results  obtained  from  the  charge  perturbation  method 
can  be  no  different  from  those  of  standard  finite-field  technique,  so  any  conclusion, 
drawn  about  the  accuracy  of  CC/MBPT  methods  will  be  the  same.  However,  the  charge 
perturbation  method  is  operationally  different.  The  charge  perturbation  method  has  been 
widely  used  by  Bishop  and  coworkers.  They  studied  simple  systems  like  LiH  (Bishop  and 
Lam,  1985),  H2  (Maroulis  and  Bishop,  1986a),  Ar  (Maroulis  and  Bishop,  1985)  and  H2O 
(Bishop  and  Pipin,  1987).  Most  of  these  calculations  have  been  performed  at  the  SCF 
level  using  extended,  carefully  prepared  basis  sets.  For  the  hydrogen  molecule  and  LiH, 
electron  correlation  was  also  included  either  by  using  an  explicitly  correlated  wavefuction 
or  the  MCSCF  method.  Both  the  SCF  and  MCSCF  methods  obey  the  Hellmann-Feynman 
theorem  and  consequently  all  first-order  moments  are  rigorously  given  as  an  expectation 
value.  Maroulis  and  Thakkar  (Maroulis  and  Thakkar,  1988,  1989)  used  the  charge 
perturbation  method  in  conjuction  with  MBPT  and  approximate  Coupled  Cluster  theories, 
to  study  polarizabilities  and  hyperpolarizabilities  of  molecules,  a topic  already  addressed 
with  standard  finite-field  method  by  Bartlett  and  coworkers  (Bartlett  and  Purvis,  1979, 
Kucharski  et  al.,  1984,  Sekino  and  Bartlett,  1986,  Pluta  et  al.,  1988).  Two  comments 
should  be  made  here.  First,  the  particular  version  of  the  Coupled  Cluster  theory  used  by 
Thakkar  namely  CCD  + ST(CCD)  or  CCD  does  not  include  single  excitations  correctly 
(or  at  all  for  CCD  !)  and  although  the  essential  orbital  relaxation  is  included  by  virtue 
of  using  CCD  with  a finite-field  Coupled  HF  reference,  the  higher-order  effects  of  single 
excitations  can  still  play  an  important  role  in  property  calculations  (Salter  et  al.,  1987) 
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and  they  are  included  in  CCSD.  More  importantly,  neither  MBPT  nor  the  usual  Coupled 
Cluster  wavefunctions  satisfy  the  Hellmann-Feynman  theorem  (Bartlett,  1987).  In  order 
not  to  introduce  errors  caused  by  the  neglect  of  the  non-Hellmann-Feynman  terms  Thakkar 
derived  polarizability  tensors  directly  from  the  perturbed  energies,  requiring  that  higher 
numerical  derivatives  be  taken  than  when  starting  from  the  dipole  moment  expansion. 
This  also  makes  the  charge  perturbation  method  loose  some  potential  advantages  for 
higher  tensors. 

Evaluation  of  First-Order  Electric  Properties 


One-electron  or  first-order  properties  can  be  evaluated  either  by  differentiating  the 
perturbed  energy  expression  (perturbation  theory  approach)  or  by  calculating  an  expec- 
tation value  of  the  electric  operator  of  interest  There  two  approaches  are  equivalent 
for  an  exact  wavefunction.  However,  for  many  approximate  wavefunctions  the  results 
are  different.  This  can  be  explained  easily  by  considering  the  general  form  of  the  first 
derivative  of  the  externally  perturbed  energy 

_ EU>  _ (*(0)|O|#(0))  .2  .... 

I dX  Ja=„-E  - (*(0)|®(0))  +s«0)P(°)-E(0)]|-sr>  (98) 

where  \I>  (0),  'I'  (A)  are  the  unperturbed  and  perturbed  wavefunction  of  the  system, 
respectively  and  O denotes  the  operator  of  interest,  while  S is  the  overlap  integral  S = 
(’3>|'I>).  The  first  term  on  the  right-hand  side  of  the  last  equation  is  the  expectation  value 
of  the  operator  O.  According  to  the  Hellmann-Feynman  theorem  (HFT)  the  remaining 
term  vanishes.  This  is  true  for  the  exact  wavefunction  or  for  completely  stationary 
wavefunctions  like  SCF  or  MCSCF  obtained  in  perturbation  independent  basis  sets. 
Neither  for  MBPT  nor  the  ordinary  CC  wavefunction  is  the  HF  theorem  satisfied. 
Alternative  UCC  methods  are  available  (Bartlett  et  al.,  1988)  that  have  this  property. 
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Let  consider  the  second  part  of  the  last  formula.  For  simplicity  neglect  the  constant 
term  2/S.  One  can  write  for  the  general  CC  wavefunction 


<*|H(0)  - E(0)|^)U=0  =(¥|H(0)  - E(0)|eT<°>|^)|>=„  + 

deTW 

(^|H(0)  — E(0)|— |$(0))|a=o 


(99) 


where  the  first  term  represents  the  derivative  of  the  SCF  reference  function.  It  is  well 
known  from  the  Coupled  Perturbed  Hartree-Fock  theory  (CPHF,  Gerratt  and  Mills,  1968) 
that  the  derivatives  of  molecular  orbitals  can  be  expressed  as  a linear  combination  of 
virtual  orbitals.  This  can  be  compactly  written  as 


|^i)|A=0  = ^x,aat.|4'(0))  (100) 

i,a 

where  x;a  is  a CPHF  expansion  coefficient  while  a*  and  i represent  now  the  creation  and 
anihilation  operators,  respectively.  The  expansion  is  valid  for  all  determinants  including 
the  ground  state.  Thus  the  last  equation  can  be  rewritten  as 

(*|[H(0)-E(0)]er|^)|A,„  = 

(*|[H(0)  - E(0)](1  + T + ir*  + ..)  £ “''WO))  = 

,,a  (101) 

£ *-<*|[H(0)  - E(0)]|D»)  + £ *(.<*•(« |[H(0)  - E(0)]|Z>“/)  + ...  = 

i, a ijab 

^c/(D,|e-THer|$(0)) 

l 

where  D denotes  the  excited  determinants  obtained  from  the  reference  function  $ by  the 
action  of  annihilation  and  creation  operators.  The  actual  form  of  the  expansion  coefficient 
c depends  on  the  level  of  excitation.  It  can  be  observed  from  the  above  equation  that  the 
only  nonvanishing  terms  on  the  rhs  of  eq.  (98)  come  from  the  excitations  neglected  in 
the  truncated  T.  This  is  due  to  the  fact  that  the  last  line  equation  is  satisfied  in  the  CC 
method.  The  hermicity  of  the  Hamiltonian  is  used  to  arrive  at  the  last  line  equation.  The 
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second  of  the  non-HFT  terms  can  be  expanded  if  we  expand  the  operator  T in  powers 
of  the  perturbation  A 


T(  A)  = T<0)  + AT*1)  + Ia2T(2)  + .... 


(102) 


However,  only  the  term  AT*1)  contributes  to  the  expansion  when  we  limit  ourselves  to 
first-order  properties.  Hence, 

deTW 

(*P(0)  - E(0)]|— ^—|*(0)>  = (¥|[H(0)  - E(0)]|T<1)£t|4.(0))  = 

E <“<1)<«,l[H(0)  - E(0)]C?)  + E i“(1)<‘<*l[H(0)  - E(0 )]Bf/)  + ...  (,03) 

i,a  ij,ab 

= J2di(Dl\e-TKeT\$) 

I 

and  the  general  formula  for  the  expansion  coefficients  can  be  written  as 


4 = E 


l7  lb 

m\n\ 


(104) 


The  expressions  for  the  coefficients  c and  d are  analogous  but  not  equivalent. 

Obviously,  in  practice  T operators  are  truncated  and  other  approximations  introduced. 
As  a consequence  of  the  truncation  at  a given  level  of  excitation  many  higher-order 
terms  remain  in  Eq.  (101)  and  (103).  The  analysis  of  these  terms  in  the  context  of 
the  CISD  method  was  discussed  thoroughly  (Nerbrant  et  al.,  1979).  The  conclusion 
from  that  paper  is  that  the  non-HFT  terms  for  CISD  contain  elements  of  the  extended 
Brillouin  matrix  between  triply  and  singly  and  triply  and  doubly  excited  determinants. 
The  connection  between  the  HFT  and  extended  Brillouin  theorem  had  been  discussed  for 
the  MCSCF  wavefunction  (Tuan,  1970).  The  use  of  a truncated  expansion  for  T does 
not  destroy  the  size-extensivity  of  the  approximate  CC  method,  however,  the  expectation 
value  ($|er  OeT|$)  becomes  an  infinite  expansion  and  further  truncation  is  necessary. 
For  the  CCSDT  model  Noga  and  Urban  (Noga  and  Urban,  1988)  retained  the  following 
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terms 

(O)  = ($1(1  + t\  + t\  + + Tl)0(  1 + Ti  + t2  + ir2T2  + T3)|$)c  (105) 

The  subscript  C indicates  that  only  connected  diagrams  are  included.  It  should  be 
emphasized  that  all  T operators  present  in  the  above  formula  are  obtained  from  fully 
converged  CCSDT-1  calculations,  therefore  they  correspond  to  a much  higher  order  of 
excitation  than  their  non-iterative  MBPT  counterparts.  From  the  iterative  solution  of  the 
CC  equations  (Kucharski  and  Bartlett,  1986)  it  follows  that  the  T2  operator  is  at  least 
first-order  in  the  correlation  perturbation,  while  with  an  SCF  starting  point  T]  and  T3  are 
at  least  of  second-order.  Sometimes  it  may  be  important  to  consider  the  contributions  to 
<0>  accurate  to  the  fourth-order  in  electron  correlation.  In  that  case  the  last  equation 
must  be  augmented  by  the  term  which  includes  the  third-order  wave  function.  This  extra 
term  is  given  as 

($|  T\t\OT2  + herm.conj.  |$)  (106) 

where  herm.conj.  stands  for  a hermitean  conjugate  of  the  previous  expression.  By  adding 
the  last  term  and  using  the  CCSDT-1  method  the  one  electron  property  is  calculated  to 
fourth-order  in  correlation  perturbation.  In  the  following  chapters  results  obtained  with 
the  second-order  wavefunction  Eq.  (105)  will  be  denoted  CCS DT- 1/2,  while  the  results 
obtained  with  the  full  fourth-order  density  will  be  denoted  CCSDT-1/4. 

Correlation  correction  to  the  one-electron  operator,  e.g.  the  dipole  moment,  are 
presented  diagramatically  in  Fig.  1.  The  diagrams  represent  the  second,  third-  and 
fourth-order  correction  to  the  first-order  property.  The  thick  blackened  line  represents  a 
T amplitude,  while  the  dotted  line  with  a black  dot  at  the  end  is  the  external  electric, 
one-electron  operator.  Diagrams  (1)  and  (2)  are  the  second-order  corrections  from  T2 
operators,  (3)  and  (4)  are  also  second-order  but  come  from  the  Ti  excitations.  Diagrams 
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Figure  1.  Diagrammatical  Representation  of  the  Electron 

Correlation  Contribution  to  the  First-Order  Property. 
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(5)  to  (8)  correspond  to  the  third-order  correction.  All  of  them  can  be  classified  as  mixed 
type:  T-JOTs,  T3OT2  and  T\OTi,  t\0T\  . Diagrams  (9)  to  (18)  represent  the  fourth- 
order  correlation  correction.  The  first  two  pairs  of  them  are  T\OT\,TlOT% , (13)  and 
(14)  are  mixed  T3OT2T2  and  hermitean  conjugates,  while  the  next  two  (15)  and  (16) 
correspond  to  the  T2T2OT2T2  terms.  The  last  two  diagrams  are  the  extra  two  terms 
added  to  make  the  full  fourth  - order  density  matrix.  Their  algebraic  expression  is  given 
in  Eq.  (106).  It  should  be  emphasized  that  the  analysis  presented  here  is  valid  only  for 
the  Hartree-Fock  reference  function. 

The  non-HFT  terms  are  difficult  to  calculate  explicitly,  especially  eq.  (103)  where 
the  first-order  perturbed  t amplitudes  are  required.  On  the  other  hand  it  is  relatively 
straightforward  to  use  a numerical  procedure  to  evaluate  these  terms.  The  difference 
between  first-order  electric  properties  calculated  as  derivatives  of  the  perturbed  energy  and 
an  expectation  value  measures  such  non-HFT  terms.  These  were  evaluated  numerically 
for  the  truncated  Cl  expansion  (Diercksen  et  al.,  1981).  The  results  for  molecules  like 
H2O,  HF  and  NH3  show  that  the  difference  in  the  dipole  polarizability  obtained  by  two 
methods  was  typically  less  than  0.01  a.u.  We  have  also  made  such  comparisons  and 
found  that  the  difference  is  sufficiently  small  for  our  purposes. 

Results  for  the  HF  molecule 


First-order  multipole  moments  calculated  with  the  [5s3p3d/3s3p]  contracted  Gaussian 
basis  set.  All  moments  reported  in  this  work  are  computed  relative  to  the  center  of  mass. 
Since,  for  linear  molecules  only  one  independent  moment  of  any  order  exists  the  irrelevant 
indices  are  omitted.  The  results  for  the  dipole,  quadrupole  and  octopole  moments  are 
given  in  Table  11.  As  can  be  seen  from  this  Table  our  correlated  results  agree  very 
well  with  the  experimental  values.  Numerical  HF  results  for  the  dipole  and  quadrupole 
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Table  11.  Electric  Moments  of  the  Hydrogen  Fluoride 


Moment 

SCF 

CCSDT-1/2 

CCSDT-1/4 

Experiment 

0.757 

0.700 

0.700 

0.707s 

0 

1.754 

1.727 

1.727 

1.72  ± 0.02b 

ft 

2.593 

2.461 

2.460 

a)  Experimental  estimate  from  Diercksen,  G.H.F.,  Kello,  V.  and  Sadlej,  A.J.  Mol.  Phys. 
49,  711  (1983) 

b)  Experimental  estimate  from  Amos,  R.D.  Chem.  Phys.  Lett.  88,  89  (1982). 


(Christiansen  and  McCullough,  1978,1979)  are  0.756  and  1.733,  respectively  again  in 
good  agreement  with  our  SCF  results.  Bishop  (Bishop  and  Maroulis,  1985)  constructed 
and  employed  8 basis  sets  ranging  in  size  from  the  contracted  [7s4p2dlf/5s2pld]  (Bl)  to 
a set  of  96  functions  [9s6p5dlf/5s4p2d]  (B8).  Their  SCF  value  for  the  dipole,  quadrupole 
and  octopole  with  the  B8  set  are  0.757133,  1.74071  and  2.62645,  respectively.  Taking 
into  account  the  difference  in  the  basis  set  size  between  B8  and  our  basis  the  overall 
agreement  at  the  SCF  level  is  very  good.  The  difference  between  CCSDT-1/2  and 
CCSDT-1/4  results  is,  as  expected,  negligible. 

The  value  of  the  octopole  moment  can  be  compared  to  the  previous  theoretical 
results.  MCSCF  pair  replacement  calculations  (Amos,  1978)  with  the  energy  optimized 
[5s3p2d/3s2p]  basis  set  yield  2.4645  a.u.,  while  the  CISD  result  (Amos,  1982)  is  2.4466 
a.u.  Both  results  are  consistent  to  our  CCSDT-1  values.  It  is  of  interest  to  present  the 
individual  correlation  contributions  to  the  values  of  the  dipole  and  quadrupole  moment. 
The  contribution  from  SCF  and  various  diagrams,  in  a.u.,  are  displayed  in  Table  12. 

As  can  be  seen  from  Table  12  the  largest  individual  contribution  comes  from  the 
mixed  terms.  These  terms  are  represented  as  diagrams  3, 4 in  the  second  order,  5 6,  7,  and 
8 in  the  third  and  13  and  14  in  the  fourth  order  of  correlation  perturbation.  Second  largest 
contribution  comes  from  the  second-order  diagrams  1 and  2.  The  additional  diagrams, 
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Table  12.  Correlation  Contributions  to  the  Dipole  and  Quadrupole  Moments 


Contribution 

0ZZ 

SCF 

0.757246 

1.753777 

<T20T2> 

-2.100642(-2) 

-1.741922(-2) 

<T2T2  OT2T2> 

2.872142  (-4) 

5.287078(-4) 

<TiT2  OT2> 

-6.535199(  4) 

- 8.782 174(-4) 

Mixed  Diagrams 

-3.2266 13(-2) 

-7.689568(-3) 

<Ti  OTi> 

-1.610971(-3) 

-2.865782(-4) 

<T3  ot3> 

-1.650825(-3) 

-1.560109(-3) 

17  and  18,  which  ensure  the  completness  of  the  fourth-order  density  matrix  contribute 
very  little.  Therefore,  one  can  expect  rather  small  differences  between  CCSDT-1/2  and 
CCS DT- 1/4  results  in  most  cases.  Both  for  dipole  and  quadrupole  the  contribution  from 
the  T2T2  is  small  and  positive. 

Table  13  displays  the  complete  set  of  the  second  - and  higher-order  tensors  as 
computed  by  the  charge  perturbation  method  from  the  correlated  perturbed  moments. 
Column  5 of  the  Table  contains  either  numerical  HF  results  or  the  SCF  results  which 
are  believed  to  be  close  to  the  HF  limit.  Most  of  the  near  Hartree-Fock  results  are  taken 
from  Bishop  (Bishop  and  Maroulis,  1985)  paper.  Only  his  B8  results  are  compared  to 
our  results.  Since  this  work  is  not  intended  to  provide  the  “best”  answers,  but  rather  to 
illustrate  the  general  accuracy  of  the  Coupled  Cluster  density  methods  we  do  not  quote  all 
SCF  nor  correlated  results  available  in  the  literature.  The  paper  of  Bishop  and  Maroulis 
can  serve  as  a rather  complete  source  of  relevant  literature  up  to  1985. 

Dipole  polarizability  a.  This  property  is  well  known  from  numerous  high-quality 
calculations.  For  parallel  components  it  is  possible  to  calculate  the  tensor  numerically 
at  the  HF  level.  The  result  obtained  by  Christiansen  and  McCullough  (Christiansen  and 
McCullough,  1979)  is  5.75  a.u.,  very  close  to  our  result  5.72.  There  is  no  numerical  HF 
result  for  the  a xx  so  we  compare  our  results  to  the  B8  SCF  value  of  Bishop.  Again,  the 
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Table  13.  Cartesian  Components  of  the  Electric  Tensors  for  HF 


Tensor 

SCF 

CCSDT-1/2 

CCSDT-1/4 

Other  SCF 

Exp.  or 
correlated 

azz 

5.72 

6.40 

6.40 

5.75a 

6.428b 

a xx 

4.49 

5.33 

5.32 

4.48  lc 

5.253b 

Pzzz 

-8.94 

-10.26 

-10.36 

-8.3a 

-9.838b 

flzxx 

0.005 

-0.54 

-0.64 

-0.4C 

-0.802b 

Az,zz 

3.98 

4.56 

4.56 

3.986c 

4.508d 

Ax.ZX 

0.76 

1.07 

1.07 

0.928c 

0.688d 

Czz,zz 

6.96 

8.16 

8.15 

7.36c 

8.63e 

Gxz,xz 

3.15 

3.71 

3.71 

4.34c 

5.15e 

Cxx.xx 

3.87 

4.70 

4.70 

5.43c 

6.53e 

Bzz,zz 

-52.6 

-64.5 

-64.5 

-56c 

Bxz,xz 

-28.2 

-40.6 

-40.5 

-35c 

Bxx.ZZ 

26.3 

41.9 

41.6 

26c 

Bxx.xx 

-39.6 

-62.4 

-62.1 

-45c 

EZ,ZZZ 

6.84 

13.30 

13.33 

6.68c 

Ex.xxx 

-0.49 

-0.11 

-0.10 

-0.76c 

7 ZZZZ 

279 

411 

413 

294c 

400b 

Txxxx 

299 

599 

597 

340c 

540b 

'Jxxzz 

120 

177 

178 

117c 

150b 

a)  Numerical  HF  results  of  Christiansen,  P.A.  and  McCullough,  E.A.,  Chem.  Phys. 
Lett.  63,  570  (1979) 

b)  CCSD  + T(4)  results  of  Sekino,  H.  and  Bartlett,  R.J.,  J.  Chem.  Phys.  84,  2726  (1986) 

c)  Near  HF  quality  SCF  results  of  Bishop,  D.M.  and  Maroulis,  G.,  J.  Chem.  Phys. 
82,  2380  (1985) 

d)  MCSCF  results  of  Amos,  R.D.,  Mol.  Phys.  35,  1765  (1978) 

e)  MBPT-SDQT(4)  results  of  Diercksen,  G.H.F.,  Kellp  and  Sadlej,  A.J.,  Mol.  Phys. 
49,  711  (1983) 


results  are  very  close.  The  CCSDT-1  results  can  be  compared  to  the  CCSD  + T(4)  results 
(Sekino  and  Bartlett,  1986)  calculated  with  an  extended  basis  set  of  [5s3p4d2f/5s3p] 
functions.  The  experimental  values  for  the  parallel  6.40  and  perpendicular  5.08  a.u. 
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polarizability  are  in  very  good  agreement  with  our  CCSDT-1  results.  The  average  value  of 
a as  given  by  CCSDT-1  is  5.68  a.u.  and  the  agreement  with  the  experimental  values:  5.60 
(Werner  and  Meyer,  1976)  and  5.52  (Diercksen  et  al.,  1983)  is  very  good.  However,  there 
is  a small  imbalance  in  the  theoretical  treatment  of  the  components  of  the  polarizability 
tensor  and  it  leads  to  an  anisotropy  of  1.08  a.u.  which  is  below  the  experimental  value 
of  1.31  ± 0.14  a.u.  (Diercksen  et  al.,  1983).  The  electron  correlation  contributes  10.6% 
to  the  total  parallel  component  and  8.1%  for  the  perpendicular  component. 

Hyperpolarizability  0.  This  third-order  property  can  be  formally  defined  as 


where  Fa  indicates  the  a Cartesian  component  of  the  external  electric  field  F.  Numerical 
HF  value  for  the  parallel  component  is  — 8.3,  which  is  above  our  — 8.94  a.u.  and  all 
SCF  results  of  Bishop  and  Maroulis.  However,  the  procedure  applied  to  calculate  the 
hyperpolarizability  was  not  numerically  very  stable,  thus  this  result  is  less  reliable  than  for 
azz.  Correlated  parallel  component  of  the  hyperpolarizability  /?  is  — 10.36  and  is  close  to 
the  — 9.838  a.u.  CCSD  + T(4)  result  of  Sekino  and  Bartlett.  The  full  fourth-order  MBPT 
value  from  the  same  paper  for  (3^  is  — 10.166.  It  should  be  noticed  that  the  inclusion 
of  triple  excitations  in  CCSD  + T(4)  is  non-iterative,  contrary  to  the  full  iterative  (but 
still  approximate  !)  scheme  of  the  CCSDT-1  method.  The  perpendicular  component 
of  the  hyperpolarizability  is  more  difficult  to  determine.  Our  SCF  calculation  yields 
the  very  small  positive  value  of  0.005  and  this  fact  illustrates  the  basis  set  deficiency, 
particulary  among  functions  perpendicular  to  the  z axis.  Among  8 SCF  results  of  Bishop 
and  Maroulis  the  range  of  values  is  very  broad:  from  — 1.7  to  +0.2  for  a quite  extended 
B7  basis  set.  CCSDT-1  results  lie  above  the  CCSD  + T(4)  and  MBPT(4)  results  of 
Sekino  and  Bartlett.  The  triple  excitation  contribution  is  reported  to  change  the  value 
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of  fixxz  from  — 0.464  to  — 0.802.  We  may  conclude  that  our  /?xxz  is  probably  slightly 
above  the  correct  value. 

Dipole-quadrupole  polarizability  A.  The  A tensor  describes  the  contribution  of  the 
field  gradient  to  the  induced  dipole  moment  and  the  contribution  of  the  field  FQ  to 
the  induced  quadrupole  moment.  The  formal  perturbation  theory  definition  of  the  tensor 
A reads 


_ yi 

nj£0 


(O[/xa|n)(n|0/?7lO)  + (Ol0flT|n)(n|/za|O) 
E„  - Eq 


(108) 


where  /i  and  0 represent  the  dipole  and  quadrupole  moment  operators.  In  this  work  the 
components  of  the  A tensor  were  computed  from  the  perturbed  quadrupole  moments.  The 
SCF  z,zz  component  3.98  is  close  to  the  best  Bishop’s  value:  3.986  a.u.  The  change  from 
the  B1  to  B8  basis  makes  the  parallel  component  smaller.  MCSCF  calculations  (Amos, 
1978)  were  performed  for  both  components  of  A.  The  agreemenet  for  parallel  components 
is  very  good.  The  electron  correlation  contribution  to  A is  relatively  small,  "12%,  and 
according  to  the  Bishop  and  Maroulis  findings  not  extremely  sensitive  to  the  quality  (of 
the  already  good)  basis  set.  For  x,zx  component  the  MCSCF  results  of  Amos  are  too  low, 
and  now  the  dynamic  correlation  becomes  a more  important  factor.  So  does  the  quality 
of  the  basis  set.  As  can  be  seen  in  Table  13  the  discrepancy  between  our  SCF  (0.76) 
and  Bishop’s  (0.928)  is  larger  than  for  the  parallel  component.  The  greater  sensitivity 
can  be  explained  by  the  fact  that  the  permanent  Oxz  is  zero,  so  we  only  deal  with  a 
rather  weak  induced  quadrupole  moment.  The  higher  angular  momentum  functions  and 
diffuse  polarization  exponents  are  needed  to  properly  describe  a situation  like  this.  Amos 
also  analyzed  vibrational  corrections  to  the  A tensor  and  predicted  a rather  substantial 
contribution  of  the  order  of  10%.  We  conclude,  the  parallel  component  is  reproduced  by 
our  CCSDT-1  calculations  rather  well,  while  the  x,zx  part  seems  to  be  too  low. 
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Quadrupole  polarizability  C.  This  quantity  is  analogous  to  the  dipole  polarizability, 
but  now  the  quadrupole  moment  is  induced  by  the  electric  field-gradient.  The  perturbation 
theory  expansion  is  now 


Cap,  y6  — yi 

n^O 


(Ol0a/?|n)(n[07g|O)  + (0lQ7gln)(n|6Q/3|0) 
E„  - E0 


(109) 


There  are  three  independent  components  of  the  C tensor  for  C0 ov  symmetry:  zz,zz  along 
the  bond,  xx.xx  perpendicular  to  the  bond  and  the  “mixed”  xz,xz  component.  Our  SCF 
values  are  consistent  with  Bishop’s  results.  The  parallel  component  shows  the  best 
agreement.  Basis  set  deficiency  leads  to  too  small  values  for  the  perpendicular  and 
“mixed”  component.  Again,  the  pertinent  permanent  quadruples  are  zero  and  we  use 
only  very  small  induced  values.  The  inclusion  of  correlation  lifts  the  parallel  component 
to  8.15  a.u.  in  good  agrement  with  the  MBPT(4)  result  (Diercksen  et  al.,  1983b).  For 
the  other  two  components  the  correlated  values  are  still  lower  than  the  SCF  results  of 
Bishop  and  Maroulis. 

Dipole-octopole  polarizability  E.  This  tensor  gives  both  the  dipole  moment  induced 
by  the  gradient  of  the  field-gradient  F^^  and  the  octopole  moment  induced  by  the 
electric  field  Fa.  The  relevant  sum-over-states  formula  is 


Ea,PiS 


E 


(0|/xQ|n)(n]»^|0)  + (0lfi/?7s|n)(n|/iQ|0) 

E„  - E0 


(HO) 


where  Q p $ denotes  the  octopole  moment  operator.  The  E tensor  may  be  of  some  interest 
for  highly  symmetrical  systems  for  which  the  octopole  Cl  is  the  first  non-vanishing  electric 
moment.  The  only  values  reported  in  the  literature  are  the  Bishop  and  Maroulis  results. 
Our  SCF  value  for  the  parallel  component  (6.84  a.u.)  is  very  close  to  their  result  (6.68 
a.u.).  The  use  of  smaller  basis  sets  leads  to  much  higher  values  of  E^zzz-  The  correlation 
contribution  is  a large  47%.  For  the  perpendicular  part  the  situation  is  less  clear.  Our 
— 0.49  a.u.  has  a smaller  magnitude  than  all  of  Bishop’s  results,  while  with  correlation, 
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one  obtains  very  small  negative  values.  Basis  set  expansion  would  likely  change  this 
component  value,  when  for  symmetry  reasons  the  permanent  moment  is  zero. 

Dipole-dipole-quadrupole  hyperpolarizability  B.  Tensor  B is  of  the  fourth-rank  and 
is  a third-order  property.  It  is  analogous  to  /?  in  a sense  that  it  has  the  same  functional 
relationship  to  the  induced  quadrupole  moment  as  (3  has  to  the  induced  dipole.  For  a 
heteronuclear  diatomic  molecule  it  has  4 independent  components.  Our  SCF  results  are  in 
general  consistent  with  the  SCF  value  of  Bishop  and  Maroulis.  For  negative  components 
our  results  are  less  negative  than  Bishop’s,  while  for  the  xx,zz  component  we  obtain  the 
same  result.  Electron  correlation  contributes:  18%,  30%,  37%  and  36%  for  the  zz,zz; 
xz,xz;  xx,zz  and  xx,xx  components  respectively.  Taking  into  account  the  difference  in 
the  basis  set  sizes  we  may  conclude  that  our  correlated  B values  are  reasonable  estimates 
of  the  experimentally  unknown  B tensor. 

Second  Hyperpolarizability  7-  This  is  a fourth-rank  tensor  and  fourth-order  property. 
The  second  hyperpolarizability  can  be  formally  defined  as 


In  principle  7 can  be  calculated  by  the  standard  finite-field  technique  like  a and  (3. 
However,  in  practice  this  method  is  extremely  difficult  to  apply  due  to  the  numerical 
problems  related  to  the  repeated  differentiation.  As  a result  the  values  of  the  second 
hyperpolarizability  7 are  sensitive  not  only  to  the  choice  of  the  basis  set,  but  particulary 
to  the  choice  of  finite-field  parameters.  For  the  parallel  component  our  SCF  results  agree 
with  the  results  of  Sekino  and  Bartlett  (280  a.u.)  and  Bishop  and  Maroulis.  The  correlated 
values  are  also  very  close  to  both  the  CCSD  + T(4)  and  MBPT(4)  values  of  Sekino  and 
Bartlett.  However,  according  to  Sekino  and  Bartlett  the  reported  results  probably  still 
represent  no  more  than  ~70%  of  the  experimental  second  hyperpolarizabilities  (Sekino 
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and  Bartlett,  1986).  The  CCSDT-1  value  of  the  perpendicular  component  yXXxx  is  higher 
than  the  analogous  result  reported  by  Sekino  and  Barlett.  Despite  the  modest  basis  set 
employed  in  these  calculations,  the  SCF  value  is  also  close  to  Bishop’s  B8  value:  340 
a.u.  and  all  the  SCF  results  of  Sekino  and  Bartlett. 

Results  for  the  H2O  Molecule 

The  electric  properties  of  the  water  molecule  have  been  intensively  studied  in  the  past. 
However,  only  recently  has  a complete  description  of  the  electric  polarizability  tensor 
been  published  (Bishop  and  Pipin,  1987).  We  used  the  charge  perturbation  method 
at  the  CCSDT-1  level  of  approximation.  The  same  basis  set  as  for  the  HF  molecule, 
[5s3p3d/3s3p]  making  a total  of  56  functions,  was  used.  The  experimental  geometry 
was  used  in  all  calculations.  Multipole  moments  were  calculated  with  respect  to  the 
center  of  mass.  The  C2  axis  coincides  with  the  z axis  and  molecule  lies  in  the  xz  plane. 
The  values  of  the  multipole  moments  are  presented  in  the  Table  14.  The  values  of 
the  zz  component  were  computed  directly  from  the  appropriate  density  matrix  and  one- 
electron  molecular  integrals.  However,  for  C2V  symmetry  only  2 independent  quadrupole 
moments  arc  possible  and  the  zz  component  can  be  related  to  the  other  components  by 

Qzz  ~ ©zz  © y y . 

The  last  relationship  reflects  the  traceless  convention  of  defining  the  multipoles 
adopted  in  this  work.  As  can  be  seen  from  Table  14  our  results  are  very  close  to  the 
experimental  data.  The  CEPA-1  calculations  with  a large  uncontracted  (Ils6p3d/5s2p) 
basis  set  (Werner  and  Meyer,  1976)  yield  0.723  a.u.  for  the  dipole  moment.  Cl  results 
for  the  xx  and  zz  components  of  the  quadrupole  moment  are  1.943  and  — 0.075  a.u. 
respectively  (John  et  al.,  1980).  An  approximate  Coupled  Cluster  method  (ACCD)  was 
used  to  calculate  many  electric  properties  of  the  water  molecule  (Dykstra,  1989).  With 
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Table  14.  Electric  Multipole  Moments  of  the  H2O  Molecule 


Moment 

SCF 

CCSDT-1/2 

CCS  DT- 1/4 

Experiment 

0.793 

0.725 

0.725 

0.7268a 

©xx 

1.921 

1.919 

1.918 

1.96  ± 0.02b 

0yy 

-1.823 

-1.837 

-1.835 

-1.86  ± 0.02b 

©zz 

-0.099 

-0.083 

-0.083 

-0.10  ± 0.02b 

a)  Clough,  S.A.,  Beers,  Y.,  Klein,  G.P.  and  Rothman,  L.S.,  J.  Chem.  Phys.  59,  2254 
(1973) 

b)  Verhoeven,  J.  and  Dymanus,  A.,  J.  Chem.  Phys.  52,  3222  (1970) 


the  contracted  [6s5p3d/3s3p]  basis  the  results  are:  — 1.880  and  — 0.056  a.u.  for  0yy  and 
0ZZ  respectively.  It  can  also  be  noticed  that  except  for  the  dipole  moment,  the  correlation 
energy  does  not  improve  the  values  of  the  moments.  As  for  the  HF  molecule  the  difference 
between  the  CCSDT-1  results  correct  through  the  second-order  in  the  wavefunction  and 
the  results  correct  through  the  full-fourth  order  in  the  corresponding  density  matrix  are 
negligible.  Table  15  displays  the  set  of  second  - and  third-order  electric  tensors  as 
calculated  by  the  charge  perturbation  method  from  the  perturbed  moments.  All  entries 
are  given  in  atomic  units. 

Dipole  polarizability  a . The  dipole  polarizability  tensor  has  been  frequently  cal- 
culated. We  confine  the  comparison  to  a few  literature  results.  Our  SCF  results  agree 
very  well  with  those  of  Werner  and  Meyer  (Werner  and  Meyer,  1976)  who  employed 
the  valence-uncontracted  basis  set  (Ils6p3d/5s2p)  with  carefully  optimized  exponents. 
A smaller,  contracted  [6s5p3d/3s3p]  basis  yields  9.077,  7.8385  and  8.3128  a.u.  for  the 
xx,  yy  and  zz  components.  MBPT-SDQ(4)  calculations  (Purvis  and  Bartlett,  1981)  with 
the  [6s5p4d/4s2p]  basis  set  give  9.87,  9.30  and  9.46  a.u.  for  the  three  components.  The 
difference  between  the  zz  and  yy  in  these  calculations  is  larger  than  in  ours:  0.16  and 
0.01  a.u.  CEPA-1  results  for  the  yy  and  zz  differ  by  0.05  a.u.  The  average  dipole  po- 
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Table  15.  Cartesian  Components  of  the  Electric  Tensors  for  H2O 


Tensor 

SCF 

CCSDT-1/2 

CCSDT-1/4 

Other  SCF 

Other 

corr./exp. 

«xx 

9.10 

10.05 

10.06 

9.24a 

9.81b,  9.87c 

Oyy 

7.91 

9.71 

9.70 

7.9  la 

9.59b,9.30c 

a ZZ 

8.42 

9.71 

9.71 

8.55a 

9.64b,  9.54c 

/^ZXX 

-9.45 

-11.2 

-11.3 

-7.87a 

-10.0C 

flzyy 

-0.39 

-4.13 

-4.41 

1.41“ 

-3.7C 

flzzz 

-6.03 

-11.2 

-11.4 

-4.68a 

-9.2C 

Ax.ZX 

3.21 

3.55 

3.55 

6.27a 

6.742d 

Ay,zy 

0.10 

0.44 

0.44 

1.40a 

1.786d 

Az,xx 

2.30 

2.31 

2.30 

1.62a 

2.54s 

Az,zz 

2.29 

2.28 

2.29 

1.68a 

2.194d, 

2.44s 

Cxx.xx 

11.6 

13.4 

13.4 

11.78a 

11.65s 

Cyy.yy 

11.5 

14.0 

14.0 

12.02a 

12.58s 

Czz,zz 

9.10 

10.6 

10.6 

10.10s 

8.76s 

Cxy.xy 

5.45 

6.91 

6.90 

8.48a 

8.48s 

Cxz,xz 

7.07 

8.28 

8.29 

11.06a 

11.33^ 

Cyz,yz 

4.47 

6.02 

5.96 

4.72a 

4.41s 

Bxx.xx 

-81.5 

-94.9 

-94.5 

-76.6  la 

Bxx.zz 

36.9 

50.2 

49.5 

39.14a 

Byy.xx 

90.3 

144.7 

144.3 

234.06s 

By,zz 

66.7 

109.8 

109.1 

67.99s 

Bzz,xx 

54.9 

82.8 

82.3 

55.78s 

BZZ,ZZ 

-94.8 

-133 

-132 

-90.91s 

BXy,xy 

-57.1 

-91.6 

-91.1 

-83.61s 

®xz,xz 

-65.0 

-88.4 

-88.2 

ed 

r-H 

00 

1 

Byz,yz 

-54.8 

-86.5 

-86.0 

-80.54s 

a)  Bishop,  D.M.  and  Pipin,  J.,  Theor.  Chim.  Acta  71,  247  (1987) 

b)  Werner,  H-J.  and  Meyer,  W.,  Mol.  Phys.  31,  855  (1976) 

c)  Purvis  III  G.D.P.  and  Bartlett,  R.J.,  Phys.  Rev.  A23,  1594  (1981) 

d)  John  I.G.,  Bacskay,  G.B.  and  Hush,  N.S.,  Chem.  Phys.  51_,  49  (1980) 
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larizability  is  experimentally  accessible.  CCSDT-1/4,  MBPT-SDQ(4),  CEPA-1  results 
are  respectively:  9.82,  9.5  and  9.68  a.u.  The  experimental  value  is  9.82  a.u.  (Werner 
and  Meyer,  1976)  perfectly  matching  our  result.  Bishop  and  Pipin  (Bishop  and  Pipin, 
1987)  used  the  MCSCF  method  (490  configurations)  with  a contracted  [5s8p3d/2s2p] 
basis.Their  values  of  the  xx,  yy  and  zz  components  are  the  following:  9.25,  9.06  and 
8.50  a.u.  It  leads  to  the  too  small  average  value  of  8.94  a.u.  Both  MBPT  and  CEPA 
calculations  were  also  performed  for  a geometry  different  from  the  experimental  one.  The 
vibrational  corrections  were  estimated  at  the  CEPA-1  level  and  were  shown  to  vary  from 
0.23  (axx)  and  0-07  (ayy)  component.  In  conclusion  we  can  say  that  the  tensor  a is  very 
well  described  at  the  CCSDT-1/4  level  with  perfect  agreement  with  the  experimental 
value  for  the  average  polarizability. 

Hyperpolarizability  0.  It  is  well  known  that  this  third-order  tensor  is  extremely 
sensitive  to  the  quality  of  the  basis  set.  The  zyy  component  was  estimated  to  vary  from 
+1.41  to — 2.76  a.u.  (Bishop  and  Pipin,  1987).  CHF  calculations  with  five  different  basis 
functions  (Lazzaretti  and  Zanasi,  1981)  give  the  range  from  1.098  to  — 0.299  a.u.  The 
value  of  — 0.5  by  Purvis  and  Bartlett  as  well  as  our  — 0.39  are  probably  too  low.  The  xxz 
component  is  reported  to  be:  — 9.5  to  — 10.6  (Lazzaretti  and  Zanasi,  1981)  for  a large 
contracted  basis  set  [9s6p3dlf/6s2pld],  — 9.6  (Purvis  and  Bartlett,  1981)  and  — 7.87  in 
the  calculations  of  Bishop  and  Pipin.  Our  value  — 9.45  fits  well  in  the  range  obtained  by 
Lazzaretti  and  Zanasi.  The  range  of  results  for  the  component  along  the  molecular  axis  is 
also  rather  broad.  Lazzaretti  and  Zanasi  report  results  as  high  as  — 3.574  a.u.,  while  the 
uncontracted  basis  of  Bishop  and  Pipin  yields  — 4.68  a.u.  Our  SCF  result  is  consistent 
with  the  previous  results.  To  conclude  the  discussion  of  the  SCF  results  it  is  interesting 
to  compare  the  results  of  analytical  Coupled  Perturbed  HF  calculations  (called  Derivative 
HF)  calculations  of  Dykstra  (Dykstra,  1989).  He  employed  a TZ  basis  augmented  with 
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polarization  functions.  The  resulting  basis  [6s5p3d/3s3p]  was  not  optimized  for  any 
particular  system  but  was  designed  to  serve  as  a medium-sized  basis  set  for  electrostatic 
calculations.  Using  this  basis  set  the  hyperpolarizability  /3  is:  — 0.5445,  — 10.029  and 
— 5.4715  a.u.  for  the  zyy,  zxx  and  zzz  components,  respectively,  in  good  agreement  with 
our  SCF  results.  The  CCSDT-1  calculations  give  values  close  to  the  previous  MBPT- 
SDQ(4)  results.  The  average  hyperpolarizability  defined  as  i(/3xxz  + f3yyz  + (3ZZZ)  is 
-15.9  for  CCSDT-1/2  and  -16.3  for  the  CCSDT-1/4.  For  comparison,  Purvis  and  Bartlett 
reported  -13.7  a.u.,  while  the  experimental  value  is  — 21.8  ± 0.9  a.u.  (Ward  and  Miller, 
1979).  Taking  into  consideration  the  notorious  difficulties  already  present  at  the  SCF 
level,  the  discrepancy  between  our  CCSDT-1  results  and  experiment  is  rather  small. 

Dipole-quadrupole  polarizability  A.  The  A tensor  was  obtained  by  the  Derivative  HF 
method  (Dykstra,  1989)  with  the  contracted  [6s5p3d/3s3p]  basis  and  by  Bishop  and  Pipin 
who  also  performed  limited  MCSCF  calculations.  The  value  of  our  x,zx  component  is 
in  fair  agreemenet  with  Dykstra’s  result  (4.0176  a.u.).  However,  our  result  is  almost  half 
of  the  corresponding  result  reported  by  Bishop  and  Pipin  (6.27  a.u.).  This  difference 
can  be  attributed  to  basis  set  limitations.  Our  basis  set  is  simply  not  flexible  enough 
to  describe  properly  the  very  small  induced  component  of  the  quadrupole  Oxz  which  is 
zero  by  symmetry  for  the  unperturbed  H2O  molecule.  The  same  effect  manifests  itself 
even  more  for  the  y,yz  component  for  which  both  the  SCF  and  the  CCSDT-1  results 
are  an  order  of  magnitude  smaller  than  the  good  quality  SCF  results  of  Bishop  and 
Pipin.  Dykstra’s  basis  yields  0.7871  a.u.,  still  far  from  Bishop  and  Pipin.  The  z,xx 
component  is  on  the  contrary  larger  than  the  values  of  Bishop  and  Pipin,  but  their  results 
suggest  that  the  better  basis  set  leads  to  smaller  values  of  Azxx.  The  DHF  result  (2.9196 
a.u)  is  probably  too  high.  The  parallel  z,zz  component  can  be  obtained  either  from  the 
perturbed  dipoles  or  perturbed  quadrupole  moments.  We  used  the  latter  formulae.  Our 
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SCF  value  lies  between  the  Bishop  and  Pipin  (1.68)  and  the  DHF  result  (2.7406  a.u.).  The 
electron  contribution  as  given  by  the  CCSDT-1  computations  is  negligible  for  the  z,xx 
and  z,zz  components  and  a moderate  “ 9.5%  for  the  Ax,zx.  For  the  small  Ay>zy  component 
correlation  changes  the  value  by  a factor  of  4.  The  MCSCF  resuts  also  show  that  for 
all  components  electron  correlation  does  not  change  the  A tensor  inordinately.  One  can 
notice  that  although  our  z,xx  and  z,zz  results  are  different  from  the  corresponding  results 
of  Bishop  and  Pipin  we  observe  the  same  approximate  equality  of  these  two  components. 
This  internal  consistency  shows  that  a big  discrepancy  for  the  y,zy  as  well  other  smaller 
differences  are  indeed  due  to  basis  set  limitations.  The  A tensor  proved  to  be  very 
difficult  to  calculate  properly.  The  few  existing  SCF  results  suggest  that  the  z,xx  and 
z,zz  components  are  probably  fairly  well  predicted.  The  situation  for  x,xz  is  worse  and 
attempts  to  calculate  the  y,yz  component  failed.  The  electron  correlation  contribution  is 
far  less  important  than  the  use  of  good  quality  basis  sets. 

Quadrupole  polarizability  C.  There  are  6 independent  components  of  this  second- 
order  tensor.  For  all  “pure”  components  our  SCF  results  are  in  very  good  agreement 
with  the  results  of  Bishop  and  Pipin.  DHF  values  are  systematically  higher:  12.8394, 
13.5368  and  11.766  a.u.  for  the  xx,xx;  yy,yy  and  zz,zz  respectively.  The  CCSDT-1 
results  are  higher  than  the  MCSCF  value  and  the  correlation  amounts  to  13%,  18% 
and  14%  for  the  xx,xx;  yy,yy  and  zz,zz  components.  The  “mixed”  components  are  all 
smaller  than  the  corresponding  results  of  Bishop  and  Pipin.  The  inflexibility  of  the  basis 
set  can  be  blamed  for  these  differences.  The  inclusion  of  electron  correlation  improves 
our  results  but  they  are  still  probably  too  small.  The  good  agreement  of  the  Cyz>y2 
component  is  fortuitous  however.  In  calculations  of  that  component  one  subtracts  the 
factorA!,)Zj,i2/2v/3  from  the  combination  of  the  perturbed  quadrupole  moments  and  our 
erroneous  small  Ay>zy  helps  to  reach  good  agreement  with  the  Bishop  and  Pipin  result. 
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The  MCSCF  correlation  contribution  is  much  smaller  than  our  CCSDT-1  results  for  this 
component.  The  DHF  calculations  yield  4.010  a.u.  for  the  yz,yz  4.776  a.u.  for  the  Cxy,xy 
and  7.014  a.u.  for  the  xz,xz  component.  All  results  are  smaller  than  our  SCF,  although 
Cxz,xz  only  slightly  and  hence  are  too  low. 

Dipole-dipole-quadrupole  hyperpolarizability  B.  To  the  best  of  the  author’s  knowl- 
edge, this  third-order  tensor  has  not  been  calculated  at  any  correlated  level.  Only  Bishop 
and  Pipin  SCF  results  can  be  compared  to  our  results.  Dykstra  uses  an  inconsistent  sign 
convention  in  his  calculations  thus  obscuring  the  comparison.  There  are  9 independent 
components  of  the  B tensor  (Dykstra  reports  12  values).  Our  SCF  results  are  consistent 
with  the  Bishop  and  Pipin  results.The  xx,xx;  zz,zz;  xx,zz;  zz,xx  and  yy,zz  components 
agree  very  well  with  the  Bishop  and  Pipin  values.  Again,  the  negative  “mixed”  compo- 
nents are  more  difficult  to  determine.  Our  SCF  values  are  higher  than  Bishop  and  Pipin 
by  approximately  30%.  The  correlation  always  lowers  the  results  beyond  the  Bishop  and 
Pipin  SCF  results.  In  conclusion  we  may  say  that  5 “pure”  components  are  probably 
accurately  predicted,  the  3 “mixed”  components  of  the  B tensor  are  fairly  well  described. 
By  extrapolating  the  correlated  results  some  reasonable  estimates  of  these  components 
can  be  made.  The  yy,xx  component  is  an  exception.  It  depends  critically  on  the  basis  set 
quality  and  our  SCF  results  amount  to  only  ~ 40%  of  the  previously  reported  value.  The 
electron  contribution  is  also  substantial  so  making  any  predictions  about  this  component 
is  not  reliable  at  this  point. 

Results  for  the  F~  Anion 

One  could  think  that  a closed-shell,  one-center  10-electron  system  like  F-  should 
not  cause  any  problem  for  modem  quantum  chemical  methods.  Despite  that  apparent 
simplicity,  however,  the  electric  properties  of  that  anion  are  extremely  difficult  to  calculate 
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correctly.  The  reason  is  the  very  diffuse  electronic  charge  distribution  demanding 
extended  basis  sets  and  the  inclusion  of  a large  portion  of  the  electron  correlation.  For 
the  charge  perturbation  method,  the  fluoride  anion  presents  a major  challenge  since  all 
calculations  must  rely  on  a proper  description  of  small  induced  moments,  permanent 
components  of  all  of  which  are  zero  for  this  spherical  system.  As  was  shown  in 
the  previous  two  sections,  the  components  for  which  the  unperturbed  corresponding 
multipole  is  zero  by  symmetry  are  usually  described  worse  than  the  other  components. 
Requirements  for  the  basis  set  made  us  use  Dykstra  [6s5p3d]  and  augmented  Dykstra 
[6s5p3dlf]  basis  sets  in  addition  to  the  standard  [5s3p3d]  set.  Since  we  are  not  interested 
in  obtaining  the  “best”  results  possible  we  decided  that  universal  (i.e.,  not  generated 
for  a specific  property  of  the  anion)  medium-sized  basis  sets  designed  for  electric 
properties  offer  a better  test  of  the  relaibility  of  the  correlated  version  of  the  charge 
perturbation  method. 

For  spherical  systems  in  the  !S  state  only  even  polarizabilities  exist:  a,  7,  C and  B. 
The  relevant  formulae  for  calculating  these  tensors  are  presented  in  the  first  section  of  the 
present  chapter.  Table  16  summarizes  our  results  and  shows  some  of  the  literature  values. 

Basis  I denotes  the  contracted  [5s3p3d]  set,  basis  II  is  the  so  called  ELP  (electric 
properties)  basis  (Liu  and  Dykstra,  1987)  [6s5p3d].  Basis  set  III  was  augmented  by  adding 
one  f function  with  exponent  0.28  optimized  previously  (Wormer  et  al.,  1982).  The  SCF 
calculations  (Maroulis  and  Bishop,  1986)  were  performed  with  two  large  basis  sets.  One 
was  a (15sl0p5dlf)  set  contracted  to  [9s6p5dlf]  (67  functions)  basis.  The  second  set  (87 
orbitals)  was  left  totally  uncontracted  and  can  be  denoted  as  (17sl0p5dlf).  The  use  of 
a larger  basis  set  usually  gives  larger  values  of  the  7 and  C tensors  and  more  negative 
values  of  B.  The  fact  reflects  the  underlying  fourth-order  Hylleraas  variational  principle, 
though  it  is  not  appllied  rigorously  to  guarantee  a bound.  Finite-field  CISD  computations 
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Table  16.  Electric  Polarizabilities  of  the  F~  Anion  in  a.u. 


Tensor, 

basis 

SCF 

CCSDT- 1/2 

CCSDT- 1/4 

Other  SCF 

Other 

corr./exp 

a basis  I 

10.21 

16.27 

16.26 

10.66a 

18.99c 

a basis  II 

10.39 

18.82 

18.82 

a basis  HI 

10.40 

18.84 

18.84 

C basis  I 

10.29 

15.09 

15.04 

17.35b 

C basis  II 

12.91 

32.23 

32.05 

23.46b 

C basis  HI 

16.82 

36.22 

36.05 

B basis  I 

-275.3 

-643.8 

-638.2 

-47  lb 

B basis  II 

-422.5 

-2157 

-2140 

-552b 

B basis  HI 

-441.2 

-2124 

-2108 

7 basis  I 

2638 

9063 

9110 

23,100d 

7 basis  II 

10,030 

103,000 

102,700 

10,900b 

79,500e 

7 basis  HI 

10,120 

99,770 

99,540 

12,900b 

a)  Numerical  CHF  result  Schmidt,  P.C.,  Weiss,  A.  and  Das,  T.P.,  Phys.  Rev.  A19, 
5525  (1979) 

b)  SCF  results  Maroulis,  G.  and  Bishop,  D.M.,  Mol.  Phys.  57,  359  (1986) 

c)  CCSDT  +T(4)  result  Kucharski,  S.A.,  Lee,  Y.S.,  Purvis  III,  G.D.P.  and  Bartlett,  R.J., 
Phys.  Rev.  A29,  1619  (1984) 

d)  CISD  result  Diercksen,  G.H.F.  and  Sadlej,  A.J.,  J.  Chem.  Phys.  76,  4293  (1982) 

e)  MBPT(4)  result  Diercksen,  G.H.F.  and  Sadlej,  A.J.,  Chem.  Phys.  131,  215  (1989) 


(Diercksen  and  Sadlej,  1982)  employed  the  contracted  [12s8p5d]  basis  derived  from 
the  (12s8p4d)  Gaussian  set  used  for  the  fluorine  atom  (Werner  and  Meyer,  1976)  by 
augmenting  it  with  several  diffuse  functions.  All  these  basis  sets  are  larger  than  ours  and 
are  individually  optimized  for  the  fluoride  anion.  Nevertheless,  as  can  be  seen  from  Table 
16  our  SCF  results  compare  well  with  the  contracted  basis  results  of  Maroulis  and  Bishop. 
However,  for  C and  B these  results  are  still  rather  far  away  from  the  best,  uncontracted 
basis  results  (“40%  for  C and  “15%  for  B).  Electron  correlation  is  very  significant  for 
all  the  tensors.  For  basis  III  it  amounts  to  45%  for  the  dipole  polarizability,  53%  for  C 
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and  almost  80%  for  the  B tensor.  Since  our  basis  sets  yield  results  which  are  far  from 
the  (unknown)  HF  limits  we  can  expect  that  the  correct  values  for  the  C and  B tensors 
may  be  rather  different  from  our  estimates. 

The  second  hyperpolarizability  is  known  to  be  very  difficult  to  determine  correctly  by 
any  quantum  chemical  method.  Numerical  differentiation  introduces  division  by  factors 
of  the  order  of  0.0054  and  thus  the  precision  of  a finite-field  calculation  is  usually  poor. 
The  charge  perturbation  method  allows  use  of  perturbed  dipoles  instead  of  the  energies 
and  is  potentially  numerically  more  stable.  Basis  I,  as  could  be  expected,  gives  very  low 
values  of  7.  Basis  II  and  III  give  results  which  are  fairly  close  to  the  expected  HF  limit 
of  '13,000  a.u.  The  correlation  contribution  is  enormous  even  by  anion  standards.  The 
CCSDT-1  values  are  almost  10  times  larger  than  the  corresponding  SCF  values.  Since 
the  charge  perturbation  method  is  not  free  from  numerical  problems  caused  by  the  choice 
of  the  charge  Q and  the  distance  R for  basis  III,  we  calculated  the  hyperpolarizability  7 
as  a function  of  Q and  R,  i.e.  of  electric  field.  The  results  are  shown  on  Fig.  2.  As  can 
be  seen  we  obtained  almost  a straight  line  and  after  extrapolation  to  the  zero  field  value, 
the  hyperpolarizability  tensor  is  still  '93,000  a.u.  For  our  largest  basis  set  [6s5p3dlf] 
we  performed  calculations  for  two  values  of  the  field.  The  results  do  not  change  much 
with  variation  of  the  field  strength  so  the  value  reported  in  Table  15  is  confirmed  for 
the  given  basis  set  and  for  the  CCSDT-1  method  of  calculation.  The  recent  MBPT(4) 
calculations  with  a large  contracted  [1  Is9p5d3f]  basis  set  yield  79,500  a.u.  for  the  second 
hyperpolarizability  (Diercksen  and  Sadlej,  1989).  It  should  be  noted  however,  that  the 
same  calculations  overestimate  the  value  of  the  dipole  polarizability  (20.58  a.u.)  by  about 
1-1.5  a.u.  It  may  be  concluded  that  the  second  hyperpolarizability  of  the  F~  anion  is 
probably  well  described  by  the  CCSDT — 1 density  matrix  method. 
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Figure  2.  Hyperpolarizability  7 of  the  F anion  in  the  electric  field 


CHAPTER  VI 


LONG-RANGE  INTERMOLECULAR  POTENTIALS 


One-Center  Multipole  Expansion 


The  intermolecular  potentials  based  on  a perturbation  theory  expansion  offer  a 
valuable  tool  to  study  systems  of  weakly  interacting  molecules.  For  many  chemically 
interesting  systems  it  is  justifiable  to  neglect  short-range  electron  exchange.  For  large 
enough  intermolecular  distances  where  in  addition  the  overlap  of  the  charge  densities 
is  very  small,  the  electrostatic  part  of  the  interaction  energy  can  be  determined  via  a 
multipole  expansion.  The  derivation  of  general  multipole  expansion  equations  is  given 
in  Chapter  II.  The  basic  energy  expression  may  be  written  in  a different  form 


The  function  E (R)  has  an  essential  singularity  at  infinity  and  henceforth  can  be  only 
asympotically  convergent.  This  means  that,  at  sufficient  large  intermolecular  distance  R, 
E (R)  can  be  arbitrarily  closely  approximated  by  a truncated  expansion 


The  Van  der  Waals  coefficients  C can  be  uniquely  determined  according  to  the  definition 


In  practice  one  does  not  know  the  function  E (R)  and  other  methods  of  determining  C 
coefficients  are  needed.  The  operator  V,  also  given  in  Chapter  II,  can  be  rewritten  to 


(112) 


k- 1 


(114) 


m= 1 
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display  the  dependence  on  R explicitly 


v = ££ 

L-0l'=0 


V , 


LV 


R(L+L'+ 1) 


(115) 


where  VLL-  contains  information  about  relative  orientation  of  molecules  and  electric 
multipole  moments  of  the  order  L and  L’.  The  net  charge  of  a molecule  is  considered 
as  a multipole  of  the  order  0.  There  are  two  different  ways  of  truncating  the  sum  in  Eq. 
(115).  The  more  natural  way  would  be  to  use  all  multipole  moments  available  for  A and 
B.  For  systems  with  known  octopoles  (L  = 3)  the  truncation  leads  to  a collection  of  terms 
with  a single  term  of  the  order  R~7.  The  alternative  expansion  is  achieved  after  collecting 
all  the  terms  with  the  same  R-dependence.  The  modified  expansion  can  be  written  as 


Vrr 


oo  A— 1 

= £E 

A=0 /=0 


Vf.A-l-1 

Rx 


(116) 


It  can  be  seen  that  for  two  neutral  molecules  interacting  through  octopoles  (L  + L + 1 = 
7)  the  last  completed  order  of  interaction  is  5.  In  other  words,  the  modified  summation 
does  not  include  some  higher  multipole-multipole  terms  but  all  orders  (in  the  distance 
R)  are  complete. 

In  the  above  formulae  the  mutipoles  of  A are  located  at  one  point  in  A,  while  the 
mutipoles  of  B are  centered  at  one  point  in  B.  The  choice  of  these  points  is  important  since 
the  mutipoles  are  not  invariant  with  respect  to  the  change  of  the  origin  of  local  coordinate 
systems.  For  example  changing  the  coordinate  origin  O to  O translated  by  a vector  r 
leads  to  the  following  expression  for  the  quadrupole  moment  in  a new  coordinate  system 


©«/?  = \ e«'(3r»«ri/?  ~ r«2<W) 


3/3  / / \ f > 1 

= ©a/?  ~ ^ocr'p  - -/i/?ra  + /x7r7<$a/?  + ~q(3rarp  - r 2Sap) 


(117) 


2 p ey  r~  fj'  a 1 r~  / -7-1*^  • ^ 

where  r'Q  is  a Cartesian  component  of  the  translation  vector.  Clearly,  for  systems  whose 
dipole  moment  is  not  equal  to  zero  the  new  quadrupole  moment  can  aquire  an  arbitrary 
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value.  Generalization  of  the  above  formula  to  the  higher  order  moments  immediately 
shows  that  only  the  first  non-vanishing  moment  is  independent  of  the  origin.  Since  the 
total  charge  is  a multipole  of  the  order  zero,  for  charged  sytems  only  the  total  electric 
charge  is  uniquely  determined. 

In  Chapter  II  the  expression  for  the  long-range  electrostatic  interactions  is  given. 
The  charge  distribution  of  the  molecule  A creates  a nonhomogeneous  electric  field, 
which  by  acting  on  the  molecule  B induces  distortion  of  its  charge  distribution.  The 
mutual  distortion  of  the  charge  distribution  is  called  induction  and  depends  upon  the 
static  molecular  polarizabilities.  For  molecules  A and  B 

AE,W  = AE^  + (118) 


and  using  the  multipole  expansion  of  the  interaction  operator  V,  Eq.  (45),  one  can 
express  the  induction  of  molecule  A as 


AEL  = - '-A*aF*Ft  - ~-C, 
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(119) 


where  FA  denotes  the  field  on  the  origin  in  molecule  A due  to  the  permanent  multipoles 
of  molecule  A.  The  electric  field  and  its  gradient  can  be  expressed  in  terms  of  the 
permanent  multipoles  of  A.  Like  in  the  case  of  electrostatic  interactions  the  formula  has 
to  be  worked  out  for  each  specific  case.  In  the  next  section  the  intermolecular  interactions 
in  the  bifluoride  ion  FHF"-  will  be  discussed.  To  illustrate  the  formalism  let  us  derive 
the  pertinent  expressions  for  the  induction.  To  simplify  the  derivation  let  assume  that  the 
anion  is  linear  and  since  the  induction  effects  are  rather  small  for  distances  of  applicability 
of  exchangeless  perturbation  theory,  only  dipole  and  quadrupole  on  HF  (A)  will  be  taken 
into  consideration.  Also,  only  the  induced  dipole  on  F-  (B)  will  be  considered,  i.e.  the 
term  linear  in  the  static  polarizability,  for  the  induced  quadrupole  moment,  proportional 
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to  the  tensor  C,  is  very  small 


(120) 


where  the  isotropic  nature  of  the  fluoride  ion  was  taken  into  account  to  drop  the  irrelevant 


suffices  for  the  polarizability.  The  electric  field  is  given  as 


ZZ^OfZZ 


(121) 


There  is  only  one  component  of  the  dipole  moment  for  the  HF  molecule  and  all  the 
off-diagonal  terms  of  the  quadrupole  moment  are  zero.  From  the  contraction  property  of 
the  T tensor  one  can  deduce  that  Taxx  + Tayy  = -Tazz  . For  heteronuclear  diatomic 
molecules  it  also  holds  that  Qxx  = Qyy  = — \QZZ  . Now  the  expression  for  the  electric 
field  on  the  fluoride  anion  reads 


Finally,  the  induction  energy  on  the  F^“  due  to  the  presence  of  the  dipole  and  quadrupole 
moments  on  the  HF  molecule  for  the  linear  orientation  is  given  as 


Now,  the  electric  field  on  A is  generated  by  a net  charge  on  B (F"~  ion).  After  some 
simple  substitution  one  arrives  at 


(123) 


The  neglected  quadrupole-induced  dipole  is  of  order  R“8,  while  the  quadrupole-induced 
quadrupole  term  is  of  order  R— 10  and  are  therefore  small.  For  the  induction  on  HF  one 
can  write  generaly 


AE 


(124) 


(125) 
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It  can  be  seen  immmediately  that  the  corresponding  induction  on  the  HF  molecule  is 
much  stronger  because  of  the  nonvanishing  charge  of  the  F-  anion.  The  usefulness  of 
the  multipole  expansion  would  be  greatly  enhanced  if  there  is  a possibility  of  improving 
the  convergence  of  the  series.  Attempts  to  do  that  are  summarized  in  the  next  section. 

Multi-Center  Multipole  Expansion 

The  main  reason  the  multipole  series  expansion  might  be  slowly  convergent  for 
long-range  and  diverges  for  shorter  distances  is  that  all  moments  on  molecule  A (B)  are 
located  at  one  point.  For  many  molecules  it  leads  to  the  necessity  of  employing  rather 
high  moments,  which  in  turn  are  not  readily  available  from  experiment  nor  from  quantum 
chemical  calculations.  Several  methods  using  a multipole  expansion  around  more  than 
one  point  have  been  proposed  recently.  One  of  them,  the  cumulative  atomic  multipole 
moments  (CAMM,  Sokalski  and  Poirier,  1983,  Sokalski  and  Sawaryn,  1987)  method 
represents  the  molecular  charge  distribution  as  a collection  of  additive  atomic  multipoles. 
CAMM  methods  are  based  on  a Mulliken  population,  at  best  a questionable  procedure, 
however,  different  definitions  of  atomic  charges  can  be  used.  An  expectation  value  of 
the  electric  operator  ak(3l •ym,  where  Greek  indices  denote  Cartesian  coordinates  and  k,l,m 
are  the  appropriate  operator  powers,  can  be  written  for  any  MO  method  as 


where  p and  q denote  atomic  orbitals  (AO).  The  first  part  is  a nuclear  contribution,  the 
second  term  is  an  electronic  contribution.  The  formula  in  the  second  line  emphasizes  the 
sum  over  atoms.  The  additive  atomic  multipoles  can  be  defined  as 


atoms 


AO  AO 
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However,  the  above  definition  suffers  from  the  dependence  of  the  multipole  upon  the 
choice  of  the  coordinate  origin.  To  ensure  an  invariance  of  the  atomic  multipoles  the 
cumulative  atomic  multipoles  were  introduced.  The  new  moments  are  obtained  from  the 
atomic  moments  by  subtracting  all  contributions  from  lower  moments.  The  first  three 

cumulative  moments  for  a particular  atom  n are  given  as 

Ha  = (a)  - qot 

©a/?  = (<*P)  - QaP  - HotP  - H0a  (128) 

= (<*/?7)  - qotP 7 - (iaPl  - HPa~l  - HlaP  - ©«/?7  ~ Q<*yP  _ © /?7a 

where  Greek  variables  denote  Cartesian  coordinates  of  a given  atom,  q is  an  atomic 

charge,  and  the  quantities  in  brackets  are  the  previously  defined  atomic  multipole  mo- 
ments. After  the  original  reference  (Sokalski  and  Poirier,  1983)  only  the  simplest  Mul- 
liken  population  analysis-based  definition  of  q is  used  in  this  work.  This  definition  can 
be  written  as 

AO 

qn  = Zn  - (129> 

pen  q 

where  the  overlap  integral  between  non-orthogonal  atomic  orbitals  is  taken  into  account. 

Adopting  this  definition  for  q the  final  expression  for  the  cumulative  atomic  dipole  is 

AO  AO 

Hat  = ((pk)®  - (pM<z»  (13°) 

p€»  q 

The  expression  for  the  other  cumulative  moments  can  be  easily  written  down.  The 
comparison  between  atomic  multipoles  and  cumulative  atomic  multipoles  for  several 
simple  molecules  like  HF,  H2O,  NH3  and  their  ions,  shows  that  the  latter  are  far  less  basis 
set  dependent,  and  less  dependent  on  small  geometry  changes.  Molecular  electrostatic 
potentials  when  computed  from  the  cumulative  moments  up  to  atomic  quadrupoles  are 
typically  about  5%  in  error  when  compared  to  the  exact  potentials  calculated  for  a SCF 
wavefunction.  On  the  other  hand,  the  potentials  calculated  from  the  atomic  multipoles  up 
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to  quadruples  can  be  almost  30%  from  a correct  value  for  distances  of  chemical  interest 
(Sokalski  and  Poirier,  1983).  The  last  characteristic  is  important  to  study  not  only  the 
long-range,  but  also  the  equilibrium  structure  of  Van  der  Waals  complexes.  The  CAMM 
method  has  been  used,  among  others,  to  model  electrostatic  effects  in  biomolecules 
(Sokalski,  1989).  The  induction  interaction,  assumed  to  be  of  an  atom-atom  type,  has 
also  been  treated  with  the  CAMM  method  (Sokalski  et  al.,  1986).  In  this  approach  the 
values  of  cumulative  multipoles  are  used  to  calculate  the  electric  field  and  its  gradient 
on  all  atoms  of  a partner.  The  CAMM  analysis  describes  the  charge  distribution  of  a 
molecule  as  a collection  of  atom  based  multipoles.  Since  the  definition  of  the  atomic 
charge  is  rather  arbitrary,  a direct  comparison  between  CAMM  values  and  experimental 
quantities  derived  from  X-ray  or  neutron  diffraction  is  not  possible. 


The  second  method  abandoning  the  traditional  one-center  multipole  model  is  the 
distributed  multipole  analysis  (DMA,  Stone,  1981,  Stone  and  Alderton,  1985).  The  basic 
idea  is  to  exploit  the  fact  that  almost  all  present  day  quantum  chemical  calculations  are 
being  performed  with  the  Gaussian  basis  sets.  The  product  of  two  primitive  Gaussian 
functions  with  exponents  a and  /?  is  itself  a new  Gaussian  located  at  the  point  vector 
P given  by 


(qA  + /3B) 
a + P 


(131) 


where  original  Gaussians  are  located  at  points  A and  B.  These  new  centers  can  be  used 
to  expand  the  molecular  charge  distribution.  It  can  be  shown  that  at  the  point  P only 
multipoles  up  to  a rank  equal  to  the  sum  of  the  angular  moment  for  the  Gaussians 
at  A and  B do  not  vanish.  In  other  words,  due  to  the  orthogonality  properties  of 
the  spherical  harmonics,  one  can  generate  a finite  number  of  new  “centers”,  and  all 
multipole  expansions  with  respect  to  these  new  centers  contain  a relatively  small  number 
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of  moments.  The  expansion  for  all  new  centers  is  an  exact  one.  For  example,  using 
only  a basis  set  of  s and  p Gaussian  functions  the  DMA  will  create  multipoles  up  to  a 
quadrupole.  Unfortunately,  the  number  of  new  points  can  grow  very  fast  due  to  the  fact 
that  all  combinations  of  two-center  products  of  primitive  Gaussians  contribute.  One  can 
reduce  the  number  of  points  by  arbitrarily  selecting  a limited  number  of  sites,  presumably 
of  some  chemical  importance  e.g.  atoms,  points  in  the  middle  of  bonds  etc.  However, 
changing  the  location  of  points  is  equivalent  to  a change  of  coordinate  origin.  The 
multipoles  contain  the  contributions  from  all  multipoles  of  lower  rank.  The  multipole 
expansion  no  longer  exactly  terminates  but  leads  instead  to  the  necessity  of  an  arbitrary 
truncation  of  the  expansion.  The  exactness  of  the  analysis  is  lost  but  the  limited  numbers 
of  points  sof  some  chemical  interest  is  usually  of  greater  value. 

It  should  be  understood  that  the  distributed  multipole  expansion  is  not  unique.  First, 
the  origin  for  the  multipole  expansion  is  arbitrary.  The  natural  choice  for  atomic  one- 
center,  i.e.  points  A and  B are  the  same,  is  the  location  of  an  atom  itself.  For  off-nuclei 
contributions  there  is  no  such  natural  choice.  Second,  multicenter  Gaussian  basis  sets 
of  modern  quantum  chemical  methods  yield  overcomplete  wavefunctions.  This  means 
that  it  is  possible  to  represent  the  basis  on  one  center  by  the  expansion  on  another. 
DMA  has  been  frequently  used  to  analyze  the  electrostatic  interactions  at  medium  and 
even  short-range  distances  for  a variety  of  the  Van  der  Waals  systems  (Buckingham  and 
Fowler,  1983). 

Both  CAMM  and  DMA  methods  use  a multicenter  expansion  of  the  molecular 
charge  distribution.  In  both  case  it  enables  extending  the  range  of  applicability  of  the 
classical,  exchangless  perturbation  theory  to  shorter  distances.  In  practical  situations  it 
improves  the  convergence  of  the  electrostatic  expansion.  Both  methods  are  approximate: 
CAMM  depends  on  the  Mulliken  definition  of  atomic  charge,  while  in  DMA  one  has  to 
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arbitrarily  select  new  origins  and  truncate  the  resulting  expansion.  The  foundations  of 
the  multicenter  multipole  expansion  were  studied  in  detail  (Vigne-Maeder  and  Claverie, 
1988)  and  a systematic  procedure  for  the  reduction  of  the  number  of  new  centers  was 
proposed.  It  was  shown  that  the  set  of  multipoles  centered  on  atoms  plus  one  point  per 
chemical  bond  exhibits  satisfactory  accuracy  for  most  cases. 


CHAPTER  VII 


APPLICATIONS 
H2O  — F^-  System 

The  mutipole  moments  and  polarizabilities  obtained  in  the  previous  chapters  can 
be  used  to  study  electrostatic  and  induction  parts  of  the  long-range  behavior  of  weakly 
interacting  complexes.  Since  the  interactions  exhibit  inherent  attractive  1/Rn  behavior 
they  are  not  suited  to  study  equilibrium  and  short-range  interactions  without  adding  a 
repulsive  component.  This  can  be  done  either  rigorously  by  using  some  version  of 
the  Symmetry  Adapted  Perturbation  Theory  as  described  briefly  in  Chapter  II,  or  using 
simpler  more  approximate  model  potential  forms. 

The  H2O  — F-  system  has  been  investigated  previously  in  the  Hartree-Fock  ap- 
proximation (Kistenmacher  et  al.,  1973).  The  basis  set  consisted  of  contracted  [7s4p2d] 
functions  for  the  F~  anion,  and  a (1  ls7pld/6s)  basis  for  H2O  contracted  to  [4s6pld/4s2p]. 
The  calculations  were  done  for  the  C2V  symmetry  of  the  system,  i.e.  the  fluoride  anion 
is  on  the  C2  axis  of  the  H2O  molecule,  and  for  a linear  configuration  where  the  F~ 
anion  is  on  a line  containing  one  of  the  two  O — H bonds  (Cs  symmetry).  The  HF 
calculations  were  of  the  supermolecule  type  and  therefore  it  was  possible  to  locate  a 
minimum.  However,  the  potential  energy  surface  was  rather  flat  around  the  minimum, 
making  results  rather  prone  to  numerical  errors.  The  interaction  energy  for  H2O  — F^“ 
was  also  calculated  using  a minimal  basis  set  using  the  CP  correction  to  reduce  the  BSSE 
(Kolos,  1979).  It  was  shown  that  the  minimal  basis  plus  the  Boys-Bemardi  correction  is 
capable  of  yielding  results  comparable  to  SCF  results  with  a larger  basis  set.  Recently 
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Table  16.  H2O  — F-  Interaction  Energy  *) 


SCF 

D(2) 

D(3) 

SDQ(4) 

SCF  a) 

SCF  b) 

4.00 

-0.81 

-4.55 

-5.44 

-4.30 

-7.22 

0.19 

4.50 

-13.74 

-15.97 

-16.81 

-15.91 

-13.25 

5.00 

-16.55 

-17.92 

-18.56 

-17.92 

-21.99 

-16.44 

5.25 

-16.43 

-17.55 

-18.07 

-17.56 

-16.51 

6.00 

-14.21 

-14.92 

-15.19 

-14.92 

-16.86 

-14.56 

7.00 

-10.82 

-11.20 

-11.30 

-11.21 

-11.24 

8.00 

-8.28 

-8.42 

-8.46 

-8.43 

-9.03 

-8.66 

10.00 

-5.22 

-5.12 

-5.13 

-5.12 

-5.52 

15.00 

-2.26 

-2.12 

-2.14 

-2.12 

-2.45 

*)  Distance  O-F  in  a.u.,  energies  in  kcal/mol 


a)  Kolos,  W.,  Theor.  Chim.  Acta  51,  219  (1979) 

b)  Kistenmacher,  H.,  Popkie,  H.  and  Clementi,  E.,  J.  Chem.  Phys.  58,  5627  (1973) 

ab  initio  calculations  at  the  CISD  level  of  theory  in  conjuction  with  an  extended  basis 
set  were  reported  (Yates  et  al.,  1988).  The  main  focus  of  the  last  paper  is  on  the  vi- 
brational properties  of  the  complex  and  therefore  a direct  comparison  with  our  results 
is  not  possible.  The  last  paper  contains  also  a rather  complete  list  of  references  to  the 
previous  work  on  this  complex. 

In  order  to  compare  the  results  of  the  long-range  multipole  expansion  we  calculated 
the  interaction  energy  by  the  supermolecule  method  using  methods  including  electron 
correlation  up  to  the  MBPT-SDQ(4).  The  BSSE  was  corrected  by  using  the  standard 
CP  technique.  Since  the  calculations  are  rather  time  consuming  a smaller  version  of  the 
basis  set  was  used  namely  [5s3p2d/3s2p].  For  selected  points  a full  basis  [5s3p3d/3s3p] 
was  employed  and  it  can  be  seen  from  the  results  the  difference  in  interaction  energy  is 
negligible.  The  results  for  the  C2V  symmetry  are  given  in  Table  16  where  the  results  of 
SCF  calculations  of  Kolos  and  Kistenmacher  et  al.  are  also  reported.  No  correction  for 
the  BSSE  effect  is  reported  in  this  Table. 
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As  can  be  seen  from  the  results  our  SCF  results  are  in  very  good  agreement  with  the 
results  of  Kistenmacher  et  al..  The  difference  around  the  energy  minimum  is  less  than  0.4 
kcal/mol.  But  for  shorter  range  the  Kistenmacher  potential  becomes  more  repulsive,  and 
at  R = 4.00  a.u.  the  difference  is  already  1.0  kcal/mol.  For  large  values  of  R our  potential 
is  about  "0.2  kcal/mol  less  attractive.  The  electron  correlation  makes  the  potential  more 
attractive  except  for  long-range  R.  The  MBPT(3)  energy  is  systematically  lower  than  the 
second  and  the  fourth-order  results.  For  large  R the  electron  correlation  becomes  less 
important,  for  R =7.00  a.u.  the  correlation  contribution  is  less  than  — 0.5  kcal/mol.  The 
results  of  Kolos  are  systematically  too  low,  at  R =4.00  uncorrected  SCF  value  is  about 
-3  kcal/mol  more  attractive  than  the  MBPT-SDQ(4).  Clearly,  the  minimum  basis  set 
yields  unacceptably  too  attractive  potential  energy  curve. 


Table  17  presents  the  H2O  — F“  potential  energy  curve  with  the  BSS  error  corrected 

by  the  standard  CP  technique.  The  results  of  Kolos  are  also  presented  for  comparison. 
Table  17.  H2O  — F-  CP-corrected  Interaction  Energy  *) 


SCF 

D(2) 

D(3) 

SDQ(4) 

SCF  a 

4.00 

-0.57 

-0.70 

-1.73 

-0.49 

-1.43 

4.50 

-13.51 

-13.25 

-14.17 

-13.19 

5.00 

-16.36 

-15.98 

-16.67 

-15.97 

-18.11 

5.25 

-16.26 

-15.87 

-16.44 

-15.86 

6.00 

-14.06 

-13.67 

-13.97 

-13.67 

-15.57 

7.00 

-10.72 

-10.36 

-10.36 

-10.36 

8.00 

-8.22 

-7.88 

-7.93 

-7.87 

-8.97 

10.00 

-5.18 

-4.90 

-4.93 

-4.89 

*)  Distance  O-F  in  a.u.,  energies  in  kcal/mol 


a)  Kolos,  W.,  Theor.  Chim.  Acta  51,  219  (1979) 


For  two  points  in  the  medium  range  the  calculations  with  a larger  basis  set 
[5s3p3d/3s3p]  were  performed.  Interaction  energies  corrected  by  the  CP  technique  are 
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Table  18.  H2O  — F~  CP-corrected  Interaction 
Energy  with  Contracted  Basis  Set  [5s3p3d/3s3p]  *) 


R=5.00 

R = 6.00 

SCF 

-16.37 

-14.09 

D(2) 

-16.05 

-13.71 

D(3) 

-16.77 

-14.04 

SDQ(4) 

-16.03 

-13.70 

*)  See  Footnotes  in  Tables  16  and  17 


displayed  in  Table  18. 

Two  conclusions  can  be  drawn  from  the  last  two  Tables.  First,  using  larger  basis  set 
does  not  change  results  much  when  the  CP  correction  is  applied.  The  SCF  energies  are 
the  same  within  the  range  to  0.03  kcal/mol.  The  largest  difference  at  the  correlated  level 
of  theory  is  — 0.10  for  the  MBPT(3)  energy.  For  R = 6.00  a.u.  the  differences  are  less 
than  0.07  kcal/mol.  The  second  conclusion  is  that  the  inclusion  of  the  intermolecular 
correlation  energy,  that  is  the  dispersion  energy,  does  not  lower  the  total  interaction 
energy.  Although,  for  the  ground  state  molecules  the  dispersion  energy  is  always  negative 
(attractive),  for  all  points  beyond  R = 4.00  a.u.  correlated  the  second  - and  fourth-order 
results  are  above  the  corresponding  SCF  values.  Third-order  MBPT  results,  however,  are 
below  the  SCF  curve  for  medium-range  points  from  4.00  to  5.25  a.u.  This  result  can  be 
explained  if  one  takes  into  account  the  differences  in  the  calculated  electric  properties  for 
H20  and  F~  at  the  SCF  and  correlated  level.  The  relevant  results  are  given  in  Chapter 
V.  The  interaction  energy  for  this  system  being  predominantly  electrostatic  is  governed 
more  by  the  changes  in  the  multipole  moments  for  the  subsystems  than  by  adding  the 
slightly  attractive  dispersion  contribution. 

Table  19  presents  the  supermolecule  calculations  for  the  Cs  geometry.  In  addition 
to  the  MBPT  calculations  the  CCSD  and  CCSD  + T(CCSD)  calculations  were  also 
performed.  Only  the  CP  corrected  values  are  presented.  As  in  the  C2v  symmetry  case 
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the  same  pattern  is  observed.  The  electron  correlation  effect  is  negligibly  small.  The  CP 
reverses  the  order  of  potential  energy  curves:  after  CP  is  applied  the  correlated  curve 
lies  above  the  SCF  curve.  The  value  of  4.68  a.u.  corresponds  to  the  minimum  reported 
by  Yates  et  al. 

The  comparison  of  the  supermolecule  results  with  the  electrostatic  and  interaction 
iterations  is  given  in  Table  20.  The  difference  between  the  approaches  decreases  with 
increasing  distance.  For  distances  equal  and  longer  than  R=6.00  a.u.  the  difference 
becomes  smaller  than  0.5  kcal/mol  at  the  correlated  level.  However,  for  R=  10.00  a.u. 
the  differences  becomes  larger  again.  For  equilibrium  distances  the  electrostatic  model 
fails  to  represent  physical  reality.  This  is  due  to  the  lack  of  any  repulsive  component 
in  the  electrostatic  model. 


Table  19.  H2O  — F-  Interaction  Energy  for  Cs 
Geometry  of  the  Complex.  Basis  Set  [5s3p2d/3s2p]*) 


4.50 

4.68 

5.50 

8.50 

SCF 

-27.47 

-27.96 

-22.97 

-7.57 

D(2) 

-26.81 

-27.59 

-22.98 

-7.39 

D(3) 

-27.60 

-28.25 

-23.24 

-7.35 

SDQ(4) 

-25.88 

-26.76 

-22.56 

-7.31 

SDQT(4) 

-25.45 

-26.45 

-22.50 

-7.31 

CCSD 

-26.62 

-27.41 

CCSD+T 

-26.35 

-27.22 

*)  Distance  O-F  in  a.u.,  energies  in  kcal/mol 


H2O  — HF  System 

The  H2O  — HF  system  is  similar  to  the  previously  studied  H2O  — F^-  complex. 
As  an  example  of  H-bonding  the  complex  has  been  frequently  studied.  Quadratic  and 
cubic  force  constants  have  been  determined  at  the  SCF  level  of  theory  using  a modified 
6-31  G basis  set  (Bouteiller  et  al.,  1980).  It  was  found  that  it  is  necessary  to  optimize 
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Table  20.  H2O  — F-  Interaction  Energy.  Supermolecule  and  Electrostatic  Interactions*) 


R 

SDQ(4)a 

SCFb 

CCSDT-l/4b 

4.50 

-13.19 

-31.85 

-29.68 

5.00 

-15.97 

-21.91 

-20.43 

5.25 

-15.86 

-19.09 

-17.91 

6.00 

-13.67 

-14.82 

-13.91 

7.00 

-10.36 

-11.28 

-10.45 

8.00 

-7.87 

-8.40 

-7.80 

10.00 

4.89 

-4.45 

-3.88 

*)  Distance  in  a.u.,  energies  in  kcal/mol 


a)  Supermolecule  calculations. 

b)  Electrostatic  and  induction  interactions 

the  geometry  of  individual  monomers  and  to  use  a good  quality  basis  set.  The  effect  of 
BSSE  has  been  investigated  (Leclercq  et  al.,  1983)  at  the  SCF  level.  It  was  shown  that  the 
minimum  of  the  C2V  arrangement  corresponds  to  — 8.13  or  — 12.59  kcal/mol  depending 
on  a basis  set  employed.  The  upper  value  is  in  good  agreement  with  an  experimental 
estimate  of — 7.8  ±1.7  kcal/mol.  The  role  of  electron  correlation  on  the  geometry  of  the 
complex  has  been  analyzed  up  to  the  MBPT(3)  level  of  theory  (Szczesniak  et  al.,  1984). 
It  was  found  that  the  overall  geometry  is  of  Cs  symmetry.  Recently,  Amos  and  coworkers 
(Amos  et  al.,  1987)  performed  the  SCF  and  MBPT(2)  calculations  on  this  complex.  The 
basis  sets  were  of  triple  zeta  quality  augmented  by  two  p functions.  The  interaction 
energy  at  the  equilibrium  is  — 5.5  kcal/mol  after  adding  zero-point  vibrational  energy. 
This  value  is  at  the  MBPT(2)  level  with  the  CP  correction  of  the  basis  set  superposition. 
The  effect  of  the  BSSE  has  been  also  been  studied  by  Latajka  and  Scheiner  (Latajka  and 
Scheiner,  1987a)  who  by  using  a modification  of  the  standard  6-31  G **  basis  showed 
that  the  BSSE  can  be  largely  reduced.  Numerous  experiments  (e.g.  Bevan  et  al.,  1980) 
have  showed  that  the  complex  H2O  — HF  is  of  Cs  symmetry  containing  a pyramidal 
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arrangement  about  the  oxygen  atom.  However,  since  the  potential  energy  is  rather  flat, 
in  this  work  only  the  C2v  symmetry  arrangement  was  treated  by  the  supermolecule-type 
correlated  approach.  Table  21  contains  the  interaction  energy  uncorrected  for  the  basis 
set  superposition.  Table  22  contains  the  interaction  energy  with  the  BSSE  reduced  by 
the  standard  CP  technique.  It  can  be  seen  from  the  results  in  Table  21  and  22  that,  like 
in  the  case  of  H2O  — F~,  the  correlation  and  the  CP  technique  reverse  the  order  of  the 
potential  energy  curves:  the  correlated  curve  lies  above  the  SCF  curve. 


Table  21.  H2O  — HF  Interaction  Energy  with  [5s3p2d/3s2p]  Basis  Set  *) 


R 

SCF 

D(2) 

D(3) 

SDQ(4) 

3.00 

-6.75 

-9.98 

-9.91 

-9.60 

3.50 

-7.70 

-10.21 

-10.23 

-10.01 

4.00 

-6.73 

-8.42 

-8.47 

-8.32 

4.50 

-5.39 

-6.51 

-6.56 

-6.45 

5.00 

-4.23 

-4.99 

-5.02 

-4.95 

5.50 

-3.33 

-3.87 

-3.89 

-3.84 

6.00 

-2.65 

-3.06 

-3.07 

-3.03 

7.00 

-1.75 

-1.98 

-1.99 

-1.96 

10.00 

-0.65 

-0.67 

-0.67 

-0.66 

*)  Distance  O-F  in  a.u.,  energies  in  kcal/mol 


Table  22.  H2O  — HF  CP-Corrected  Interaction  Energy. 


SCF 

D(2) 

D(3) 

SDQ(4) 

3.50 

-7.44 

-7.38 

-7.43 

-7.23 

4.00 

-6.49 

-6.35 

-6.41 

-6.28 

4.50 

-5.22 

-5.03 

-5.08 

-4.99 

5.00 

-4.10 

-3.90 

-3.93 

-3.87 

6.00 

-2.56 

-2.38 

-2.40 

-2.36 

7.00 

-1.71 

-1.59 

-1.61 

-1.59 
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The  difference  in  SCF  values  of  the  electrostatic  moments  and  polarizabilities  and  the 
correlated  values  of  the  moments  provides  an  explanation.  From  the  ab  initio  calculations 
is  known  that  the  H2O  molecule  is  slightly  deformed  in  the  complex  and  that  the  bond 
length  in  the  HF  molecule  is  different  from  that  in  the  free  molecule.  In  order  to  obtain 
very  good  agreement  with  an  experimental  vibrational  data  it  is  necessary  to  fully  optimize 
all  bond  lengths  and  angles  for  all  O — F distances.  In  order  to  assess  an  error  made 
by  using  [5s3p2d/3s2p]  basis  set  the  calculations  were  performed  for  five  points  using  a 
larger  Gaussian  basis  set  [5s3p3d/3s3p].  The  differences  in  energies  calculated  without 
the  BSSE  correction  are  smaller  than  0.25  kcal/mol.  The  CCSD  and  CCSD  + T(CCSD) 
calculations  were  made  for  R = 3.50  a.u.  The  results,  after  the  CP  corrections  are:  — 7.26 
and  — 7.33  kcal/mol,  respectively,  to  be  compared  to  — 7.23  obtained  at  the  SDQ(4)  level. 
The  results  indicate  that  the  basis  set  related  errors  are  likely  to  be  more  important  than 
the  inclusion  of  some  excitations  to  an  infinite  order.  Finally,  instead  of  the  standard  CP 
technique  a simple  modfication  was  used.  Only  the  diffuse  and  polarization  functions 
located  on  a neighbor  were  taken  as  a basis  set  for  the  calculation.  The  calculations 
were  done  for  one  point  R = 3.50  a.u.  and  the  results  are  given  in  Table  23.  The  last 
column  entries  are  differences  between  the  standard,  or  full  CP  and  the  approximate  CP 
corrections.  It  can  be  seen  that  an  approximate  virtual  only  version  of  the  CP  correction 
is  smaller  by  0.5  — 0.6  kcal/mol.  It  amounts  to  about  15%  of  the  standard  CP  correction. 
However,  as  can  be  seen  from  the  Table  23  this  small  difference  can  change  qualitatively 
the  description  of  the  interaction.  The  fractional  CP  leads  to  the  more  attractive  correlated 
curve,  while  the  full  CP  corrections  change  the  order  of  the  curves.  For  the  short  distance 
of  R = 3.50  a.u.  the  dispersion  effects  are  negligible  so  the  difference  is  caused  entirely 
by  the  description  of  the  charge  distribution  and  exchange  effects.  Both  these  effects 
are  better  described  in  the  full  dimer  basis,  hence  the  full  CP  results  are  probably  more 
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Table  23.  H2O  — HF  Interaction  Energy.  Standard  and  Approximate  CP  Corrections. 


Method 

No  CP 

Full  CP 

Approx.  CP 

Difference 

SCF 

-7.77 

-7.44 

-7.49 

0.05 

D(2) 

-10.21 

-7.38 

-7.94 

0.56 

D(3) 

-10.23 

-7.43 

-7.98 

0.55 

SDQ(4) 

-10.01 

-7.23 

-7.76 

0.53 

SDQT(4) 

-10.19 

-7.27 

-7.85 

0.58 

CCSD 

-10.04 

-7.26 

-7.80 

0.54 

CCSD  + T 

-10.25 

-7.33 

-7.91 

0.58 

relaible. 


FHF^-  System 


For  the  hydrogen  bifluoride  FHF-  system  only  the  multipole  analysis  was  performed. 
The  electronic  structure  and  vibrational  spectrum  of  this  anion  has  been  evaluated  recently 
by  the  CCSDT-1  method  using  a contracted  [5s3p2d/3s2p]  basis  set  (Spirko  et  al.,  1989). 
These  calculations  were  of  supermolecule  type  and  were  focused  on  the  near-equilibrium 
region  of  the  potential  energy  surface.  Previous  CISD  calculations  reported  rather  large 
changes  of  the  vibrational  frequencies  with  the  level  of  theory  and  the  basis  set  used 
(Janssen  et  al.,  1986).  The  geometry  of  the  complex  and  vibrational  frequencies  have 
been  studied  at  the  MBPT(2)  level  of  approximation  employing  various  basis  sets  of  6-31 
quality  (Frisch  et  al.,  1986).  The  bonding  energy  was  obtained  without  any  correction  for 
the  basis  set  superposition  error.  Potential  and  dipole  moment  surfaces  and  corresponding 
vibrational  frequencies  and  transition  dipole  moments  were  obtained  from  the  CEPA-1 
calculations  (Botschwina,  1987).  From  the  early  ab  initio  SCF  calculations  (Almlof, 
1972)  it  is  known  that  the  equilibrium  structure  is  linear,  so  it  is  appropriate  to  start  from 
this  geometry. 

The  fluoride  anion  does  not  possess  any  non-vanishing  electric  multipole  moments  except 
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for  the  total  charge.  It  makes  the  electrostatic  interaction  expression  particulary  simple 

AEe/ec<r  = Q(jp  + ^3  + #)  (132) 

where  Q denotes  the  total  charge  of  the  fluoride  ion  ( — 1),  and  for  a linear  geometry  of 
the  complex  //,  0 and  O should  be  understood  as  the  appropriate  diagonal  components. 
The  expressions  for  induction  are  given  in  Chapter  VI  and  need  not  to  be  repeated  here. 
Table  23  gives  the  total  electrostatic  contribution  as  a function  of  R for  the  linear  FHF- 
anion.  As  can  be  seen  from  the  data  the  presence  of  a noncompen sated  charge  of  the 
fluoride  ion  causes  very  strong  attractive  interaction.  Also,  as  a direct  consequence  of 
smaller  values  of  the  correlated  multipoles  than  the  corresponding  SCF  values,  the  SCF 
energies  are  stronger  from  ~ 2.9  to  0.4  kcal/mol.  One  can  also  notice  the  change  in 
the  composition  of  electrostatic  energy:  at  short  distances  the  dipole-charge  interaction 


amounts  to  "55%  of  the  interaction,  for  R = 7.00  the  contribution  is  already  70%. 
Table  24.  Electrostatic  Energy  of  FHF^*) 


R in  a.u. 

SCF 

CCSDT-1/4 

4.00 

-53.27 

-50.44 

4.25 

-45.65 

-43.22 

4.50 

-39.52 

-37.37 

4.75 

-34.54 

-32.63 

5.00 

-30.43 

-28.73 

5.50 

-24.11 

-22.74 

6.00 

-19.56 

-18.42 

7.00 

-13.59 

-12.78 

8.00 

-9.98 

-9.36 

10.00 

-6.02 

-5.63 

*)  Distances  from  the  center  of  mass  in  a.u.,  energies  in  kcal/mol 


Table  25  contains  the  induction  contribution  for  the  linear  FHF^- complex.  Total 
induction  is  the  sum  of  the  induction  on  the  F^“  anion  and  induction  on  the  HF  molecule. 
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The  distance  R is  measured  from  the  center  of  mass  of  the  hydrogen  fluoride.  All 
energies  are  in  kcal/mol.  As  can  be  seen  from  the  Table  induction  is  much  weaker  than 
the  electrostatic  energy.  Nevertheless,  it  cannot  be  ignored.  For  R = 4.50  a.u.  induction 
is  about  29.5%  of  the  electrostatic  interaction,  even  for  R = 8.00  the  induction  is  still 
6.5%  of  the  electrostatic  interaction.  The  induction  on  the  HF  molecule  is  stronger  than 
the  induction  on  the  F^“  anion.  This  was  already  evident  from  the  R-dependence  of  the 


expressions  in  Chapter  VI. 

Table  25.  Induction  Energy  of  FHF^-  *) 


R in  a.u. 

SCF  a 

CCSDT-l/4b 

CCSDT-l/4c 

CCSDT-l/4d 

4.00 

-10.42 

-11.80 

-9.05 

-20.85 

4.25 

-7.98 

-9.025 

-5.94 

-14.965 

4.50 

-6.21 

-7.02 

-4.02 

-11.04 

4.75 

-4.905 

-5.54 

-2.76 

-8.30 

5.00 

-3.93 

-4.43 

-1.94 

-6.37 

5.50 

-2.60 

-2.935 

-1.01 

-3.945 

6.00 

-1.79 

-2.02 

-0.56 

-2.58 

7.00 

-0.93 

-1.05 

-0.20 

-1.25 

8.00 

-0.53 

-0.60 

-0.08 

-0.68 

10.00 

-0.21 

-0.23 

-0.02 

-0.25 

*)  Distances  in  a.u.,  energies  in  kcal/mol 


a)  Induction  on  the  HF  molecule 

b)  CCSDT-1/4  induction  on  the  HF  molecule 

c)  CCSDT-1/4  induction  on  the  F-  anion 

d)  Total  CCSDT-1/4  induction 

The  angular  dependence  of  the  electrostatic  and  induction  interactions  can  be  worked 
out  from  the  angular  dependence  of  the  T tensors.  General  forms  of  the  electrostatic 
interaction  formulae  have  been  tabulated  for  arbitrary  geometries  of  the  systems  (Price 
et  al.,  1984).  For  induction  and  approximate  dispersion  the  angular  dependence  has  been 
analyzed  in  terms  of  irreducible  tensors  (Leavitt,  1980).  In  our  case  the  final  formulae 
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are  quite  simple 


AEejec<r  — Q 


acosO  0(3 cos2 6 - 1)  ft(5co.s30  - 3 cosO) 

~W~+  2B?  + 2R4 


(133) 


for  the  electrostatic  energy  and 


(134) 


for  the  induction.  An  angle  9 is  between  the  bond  in  HF  and  the  straight  line  from 
the  center  of  mass  to  the  F~  ion,  a is  the  average  of  the  dipole  polarizability  tensor 
and  so  is  C.  The  rest  of  the  symbols  have  their  usual  meaning.  The  presence  of  the 
Legendre  polynomials  in  the  electrostatic  formula  can  be  readily  recognized.  The  angular 
dependence  is  illustrated  in  Fig.  3.  where  the  total  induction  is  plotted  as  a function 
of  9 for  a constant  R equal  to  5.00  a.u.  Fig.  4 displays  the  angular  dependence  of  the 
electrostatic  interactions.  Three  R distances  are  shown:  5,  6 and  7.00  a.u.  Also,  for  a 
sake  of  comparison  the  total  induction  energy  curve  for  R = 5.0  a.u.  is  displayed.  First  of 
all,  it  is  obvious  that  the  angular  dependence  of  electrostatic  contribution  is  much  stronger 
than  the  analogous  induction  dependence.  But,  there  is  an  angle  for  which  fast  changing 
electrostatic  contribution  is  equal  to  slowly  varying  induction.  Also,  for  the  angles  around 
60  degrees  the  induction  is  the  only  attractive  force  for  the  FHF-  complex.  The  second 
conclusion  from  Fig.  4 is  that  the  angular  dependence  of  the  electrostatic  interaction  is 
very  weak  for  longer  R.  Strong  angular  dependence  of  the  electrostatic  interactions  can 
be  used  to  determine  the  shape  of  H-bonded  sytems  in  equilibrium. 

Standard  multipole  expansion  cannot  be  applied  when  the  distance  between  two 
systems  becomes  too  short.  In  order  to  enhance  the  applicability  of  classical  electrostatic 
ideas  the  CAMM  method  was  used.  As  can  be  seen  from  the  Table  26  the  difference 
between  the  CAMM  and  standard  analysis  becomes  evident  only  at  short  distances.  The 
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Figure  3.  Induction  Energy  as  a Function  of  6 for  the  FHF~  System 
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Figure  4.  Electrostatic  and  Induction  Energy  as  a Function  of  6 for  the  FHF-  System 
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Table  26.  Electrostatic  Interaction  of  the  FHF-  System  using  CAMM  Method 


R in  a.u. 

SCFa 

Ab 

CCSDT-1C 

Ad 

4.00 

-41.62 

11.65 

-39.29 

11.15 

4.25 

-37.54 

8.11 

-35.29 

7.93 

4.50 

-34.27 

5.25 

-32.36 

5.01 

4.75 

-31.72 

2.82 

-30.07 

2.56 

5.00 

-29.41 

1.02 

-27.80 

0.93 

5.50 

-23.91 

0.20 

-22.53 

0.21 

6.00 

-19.58 

0.02 

-18.40 

0.02 

a)  SCF  using  Cumulative  Atomic  Multipole  Moments  Method 

b)  The  difference  between  standard  multipole  and  CAMM  method 

c)  CCSDT-1  result  using  CAMM  method 

d)  The  difference  between  CCSDT-1  and  CCSDT-1  CAMM  method 


electrostatic  part  of  interaction  energy  can  be  calculated  easily  once  the  CAMM  are 
known.  However,  induction  is  more  complicted  to  determine.  In  this  work  the  simple 
approach  was  adopted:  atomic  multipoles  are  used  as  sources  of  electric  field  affecting 
the  multipoles  of  a neighbor.  This  approach  introduces  an  atom-atom  potential  method 
widely  used  to  study  properties  of  molecular  crystals,  large  organic  molecules  etc.  (e.g. 
Pertsin  and  Kitaigorodsky,  1987). 

For  distances  longer  than  5.00  a.u.  the  standard  multipole  analysis  and  the  CAMM 
procedure  give  results  differing  by  less  than  1.0  kcal/mol.  However,  for  shorter  distances 
the  CAMM  results  are  significantly  less  attractive  than  the  standard  expansion.  Since  the 
attractiveness  of  electrostatic  interactions  at  short  distances  sharply  contradicts  physical 
reality  one  may  conclude  that  the  CAMM  provides  much  more  reasonable  description  of 
intermolecular  forces.  For  systems  dominated  by  the  electrostatic  interactions  the  CAMM 
method  enhances  the  range  of  applicability  of  electrostatic  models  by  25%.  However,  the 
inclusion  of  induction  is  less  rigorous:  by  the  use  of  atom-atom  potentials.  For  strongly 
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interacting  hydrogen-bonded  systems  the  gain  in  better  description  of  the  electrostatic 
interactions  recompensates  for  this  drawback. 


CHAPTER  VIII 


CONCLUSIONS 

The  multipole  moments  and  related  polarizabilities  for  selected  simple  molecules  have 
been  calculated  using  a density  matrix  Coupled  Cluster  method.  This  approach  permits 
the  efficient  and  accurate  computation  of  electric  moments  as  the  expectation  value  of 
relevant  one-electron  operators.  At  the  highest  level,  Coupled  Cluster  method  contains 
converged  singles  and  doubles  and  the  most  important  triple  excitation.  The  inclusion 
of  Ti  is  sufficient  to  introduce  most  orbital  relaxation  effects.  It  is  therefore  shown 
that  the  non-Hellmann-Feynman  terms  expressed  as  higher-order  generalized  Brillouin 
matrix  elements  are  likely  to  be  small.  The  overall  agreement  with  experiment  and  other 
high-level  quantum  chemical  calculations  is  very  good.  The  availability  of  an  accurate 
CC  dipole  and  other  moments  permits  evaluating  the  components  of  electric  tensors  by 
differentiation.  Several  test  calculations  indicate  that  the  differentiation  of  the  induced 
dipole  moment  is  numerically  more  stable  than  the  ordinary  finite-field  technique.  For  a 
majority  of  electric  properties  the  difference  between  the  results  obtained  with  the  second- 
order  wavefunction  and  fourth-order  density  is  shown  to  be  negligible.  The  advantage 
of  the  CC  density  matrix  approach  in  conjuction  with  the  charge  perturbation  is  that  its 
applicability  can  be  rather  easily  extended  by  adding  higher  multipoles  to  include  new 
higher-order  properties.  In  general,  the  components  of  electric  tensors  which  can  be 
related  to  the  non-vanishing  component  of  the  permanent  multipole,  e.g.  z component  of 
the  dipole  moment  of  HF,  are  better  described  than  the  components  related  to  a vanishing 
component  of  a multipole,  e.g.  all  moments  of  the  Ne  atom.  The  modest  size  of  the  basis 
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set  is  responsible  for  this  behavior.  For  the  flouride  anion  our  calculations  yield  much 
larger  values  for  the  second  hyperpolarizability  tensor  than  those  previously  reported  by 
Cl,  finite-field  type  calculations.  In  this  case  all  tensors:  a,  C,  B and  7 are  related  to 
vanishing  multipoles  of  F- . 

It  is  shown  that  the  basis  set  superposition  error  which  is  known  to  change  dramat- 
ically the  interaction  energy  is  of  much  less  importance  for  the  dipole  polarizability  of 
interacting  Ne  atoms.  The  calculation  of  the  polarizability  function  of  Ne2  displays  a 
shallow  local  minimum  in  the  anisotropy  curve.  A similar  effect  has  not  been  observed 
so  far  for  the  other  noble  gas  dimers. 

The  multipole  moments  and  polarizabilities  have  been  used  to  construct  several 
intermolecular  potentials  for  hydrogen-bonded  species.  For  such  cases  the  electrostatic 
forces  dominate  the  interaction,  therefore  the  omission  of  the  dispersion  does  not  cause 
a much  error.  The  long-range  behavior  for  the  H20  — F system  is  compared  with  the 
supermolecule  calculations.  For  the  FHF“  anion  only  the  electrostatic  plus  induction 
models  are  considered.  The  induction  is  usually  much  weaker  for  hydrogen-bonded 
sytems,  nevertheless,  it  becomes  important  at  long  range.  The  Cumulative  Atomic 
Multipole  Moment  method  has  been  used  to  improve  the  performance  of  the  classical  or 
standard  multipole  expansion.  The  numerical  differences  between  standard  and  CAMM 
potentials  begin  to  appear  at  distances  close  to  the  equilibrium  distance.The  differences 
between  SCF  and  highly  correlated  electrostatic  potentials  are  found  to  be  rather  small.  It 
can  be  also  seen  from  the  correlated  results  for  the  H20  — F-  and  F120  — HF  systems 
that  the  attractive  dispersion  part  is  smaller  than  the  difference  between  the  electrostatic 
models  at  the  SCF  and  the  correlated  level. 
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